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The  complex  variable  boundary  element  method  (CVBEM)  for  simply- 
connected  domains  is  extended  in  three  new  areas.  A  complete  formulation  of 
the  CVBEM  using  quadratic  elements  is  presented.  The  derivation  follows  the 
format  for  linear  elements  given  in  the  literature,  and  both  the  nodal-  and 
interior- point  equations  axe  given  in  detail.  Second,  a  methodology  for  the 
treatment  of  corners  in  the  CVBEM  is  developed  for  all  possible  combinations  of 
Dirichlet,  Neumann,  and  Robin  boundary  conditions  at  the  corner.  It  is  shown 
that  no  special  augmentation  of  the  CVBEM  matrix  equations  is  necessary  when 
one  follows  the  guidelines  set  forth  in  this  work.  Finally,  an  implicit,  iterative 
variation  of  the  linear-element  CVBEM  is  applied  to  nimierical  grid  generation 
in  two-dimensional,  simply-connected  spatial  domains. 

The  quadratic-element  formulation  is  successfully  tested  by  solving 
example  problems  with  available  analytical  solutions.  Analytical  solutions  are 
also  used  to  verify  the  corner  handling  methodology,  and  the  CVBEM  results  axe 
compared    with    those    from    a    real    variable    boundary    element    method    for 

vi 


accuracy.  The  application  of  the  CVBEM  to  numerical  grid  generation  is 
demonstrated  by  generating  grids  for  several  irregular  geometries.  The  quality  of 
the  resulting  grids  is  assessed  by  visual  inspection  and  by  examining  the 
distribution  of  the  metrics. 
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CHAPTER  I 
INTRODUCTION  AND  LITERATURE  REVIEW 


Over  the  past  ten  years,  a  new  class  of  numerical  methods  for  solving 
partial  differential  equations  (PDEs)  has  risen  to  challenge  the  well-established 
finite-difference  methods  (FDMs)  and  finite-element  methods  (FEMs)  [1,2]. 
Unlike  the  FDMs  and  FEMs  which  require  that  the  entire  domain  be  broken 
down  into  grid  points  or  elements  (discretized),  these  relatively  new  methods, 
known  as  boundary  element  methods  (BEMs),  only  involve  discretization  of  the 
boundaries  of  the  domain.  In  this  way,  the  BEMs  effectively  reduce  the 
dimension  of  the  discretization  by  one,  and  in  the  process,  significantly  decrease 
the  number  of  simultaneous  equations  which  need  to  be  solved.  These  benefits 
do  not  come  free,  however,  as  the  matrices  associated  with  these  systems  of 
equations  tend  to  be  neaxly  fully  populated  in  the  BEMs,  while  they  axe 
comparably  sparse  in  the  FEMs  and  the  implicit  FDMs.  Even  with  this  tradeoff, 
it  has  been  shown  that  the  BEMs  are  often  much  more  computationally  efficient 
than  the  FDMs  or  FEMs  [1-6]. 

One  way  of  classifying  the  different  kinds  of  BEMs  is  according  to  the 
variables  used  in  their  formulation.  Traditionally,  the  BEMs  have  been 
formulated  using  real  variables  (RVBEMs),  but  a  recent  and  powerful  advance 
involves  the  use  of  complex  variables,  the  result  being  a  method  known  as  the 
complex  variable  boundary  element  method  (CVBEM)  [7].  This  new  method  is 
limited  in  application  to  Laplace's  and  Poisson's  equations,  and  is  further  limited 
to    two-dimensional    (2-d)    domains,    but    it    does    possess    the    following    two 


significant  advantages  over  the  RVBEMs:  1)  approximations  are  made  only  at 
the  boundary  of  the  domain,  and  2)  all  integrations  are  carried  out  analytically 
without  the  need  for  numerical  integration.  These  two  factors  combine  to  give 
the  CVBEM  excellent  potential  for  both  accuracy  and  efficiency. 

The  origins  of  the  CVBEM  can  be  traced  to  a  paper  by  Himt  and  Isaacs  [8] 
where  it  was  applied  to  the  solution  of  groundwater  flow  problems.  Hromadka 
and  Guymon  [9]  subsequently  used  the  method  to  predict  freezing  fronts  in  soils 
and  then  formalized  the  rather  loose  development  into  the  CVBEM  [10].  The 
method  was  also  shown  by  Hromadka  [11]  to  be  a  generalization  of  the  analytic 
function  method  of  Van  Der  Veer  [12].  In  the  process  of  refining  the  CVBEM, 
Hromadka  [13]  developed  a  technique  for  error  visualization  at  the  domain 
boundary  and  determined  relative  error  bounds  for  the  method  [14].  He  also 
considered  the  use  of  variable  trial  functions  for  improving  accuracy  [15]  and 
investigated  the  proper  placement  of  collocation  points  on  the  domain  boundary 

[16].  : 

Various  physical  phenomena  have  been  modeled  by  the  CVBEM  including 
advective  and  groundwater  contaminant  transfer  [8,17,18],  conduction  heat 
transfer  [16],  prediction  of  freezing  fronts  in  soil  [9,19],  and  stratified  flows  [20]. 
All  of  these  applications  and  developments  have  been  detailed  in  a  book  on  the 
CVBEM  by  Hromadka  and  Lai  [7];  however,  more  recent  developments  in  the 
CVBEM  (post  1987)  are  not  covered  in  this  book.  The  CVBEM  has  since  been 
extended  from  simply-connected  to  doubly-  and  multiply- connected  domains  by 
Kassab  [21],  Kassab  and  Hsieh  [22],  and  Hsieh  and  Kassab  [23].  It  was  found 
that  the  complex  potential  along  cuts  in  the  domain  does  not  cancel  out  but 
results  in  a  complex  stream  function  that  plays  the  role  of  a  perturbation  in  the 
nodal  equations.  About  the  same  time,  Harryman  et  al.  [24]  and  Hromadka  [25] 
applied    the    CVBEM    to    specific    multiply-connected    domains    where    such 


perturbations  did  not  appear  due  to  a  zero-flux  condition  on  the  interior 
boundaries.  Very  recently,  Mokry  [26]  has  applied  the  CVBEM  to  external 
potential  flows.  Examples  were  given  for  flows  over  airfoils  whose  exact  solutions 
are  known.  In  all  cases,  the  agreement  of  the  CVBEM  results  with  the  theory 
was  excellent. 

This  dissertation  is  concerned  with  extending  the  CVBEM  in  three  areas, 
namely,  boundary  corner  treatment,  higher-order  elements,  and  numerical  grid 
generation.    Each  of  these  topics  will  be  discussed  in  the  next  three  sections. 

1.1  The  Corner  Problem 

An  interesting  problem  common  to  many  BEMs  is  the  proper  treatment  of 
the  corners  or  edges  in  the  boundaries  of  the  domain.  A  corner  (in  2-d)  or  edge 
(in  3-d)  is  defined  as  a  point  or  locus  of  points  where  the  tangent  to  the 
boimdary  possesses  a  sharp  discontinuity.  This  discontinuity  must  be  due  solely 
to  the  problem  geometry  and  not  due  to  any  discretization  scheme  associated 
with  the  BEM.  At  such  a  corner  or  edge,  quantities  associated  with  the 
direction  normal  to  the  boundary  are  double- valued,  and  the  conventional 
RVBEMs  can  produce  systems  of  equations  which  are  underconstrained. 

In  2-d  potential  problems,  the  three  traditional  boundary  conditions  are  1) 

the  Dirichlet  condition  where  (j)  is  known  and  o-  is  unknown,  2)  the  Neumann 

condition  where  -^  is  known  and  (f)  is  unknown,  and  3)  the  Robin  (or  convective) 
on 

condition  where  both  6  and  -^  are  unknown.     The  symbol  (f)  represents  the 

on 

potential  function,  while  n  represents  the  direction  normal  to  the  boundary.  The 
application  of  these  three  conditions  logically  leads  to  nine  possible  boundary- 
condition  combinations  at  any  corner;  see  Table  1-1.  Here,  the  subscripts  U  and 
D  refer  to  the  "upstream"  and  "downstream"  sides  of  the  corner,  respectively, 
with  the  direction  of  travel  along  the  boundary  always  placing  the  domain  to  the 
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TABLE    1-1.      The   nine    traditional   boundaxy   condition   combinations 
at  a  corner  node  in  a  potential  problem. 


Case 


Upstream  Nodes  Downstream  Nodes 

Condition    Quantity  Specified     Condition    Quantity  Specified 


1 

Dirichlet 

<l> 

Neumann 

d<i> 

2 

Neumann 

Dirichlet 

<i> 

3 

Neumann 

d<f> 

Neumann 

4 

Dirichlet 

^ 

Robin 

none 

5 

Rxabin 

none 

Dirichlet 

<!> 

6 

Robin 

none 

Robin 

none 

7 

Neumann 

Robin 

none 

8 

Robin 

none 

Neumann 

9 

Dirichlet 

<i> 

Dirichlet 

<!> 

left.  Also,  the  symbol  s  will  be  used  to  represent  the  direction  tangent  to  the 
boundary.  For  the  optimum  situation  where  both  the  upstream  and  downstream 
conditions  axe  available  at  the  corner  node,  the  unknowns  and  corresponding 
number  of  equations  which  can  be  applied  in  ordinary  RVBEM  formulations  are 
given  in  Table  1-2.  In  all  nine  cases,  one  boimdary  element  nodal  equation  is 
always  available  for  application.  Thus,  with  only  one  unknown.  Cases  1  through 
3  pose  no  problem.  Cases  4  and  5  have  two  unknowns,  but  in  both  cases  an 
extra  equation  is  available  from  the  Robin  condition.  Case  6  has  three 
unknowns,  but  two  additional  equations  are  available  from  the  Robin  conditions. 
Cases  7  and  8  have  two  unknowns  each,  but  once  again,  the  Robin  condition  can 
be  used  as  an  extra  equation  in  either  case.  Finally,  Case  9  has  two  unknowns, 
but  no  additional  equation  is  available  for  application.  This  renders  Case  9 
underconstrained,  and  it  is  here  that  investigators  have  been  forced  to  seek 
methods  for  special  treatment  of  corners  in  the  RVBEMs. 

Probably  the  simplest  way  to  treat  Ccise  9  is  to  assume  that  (^1/  and 
i-^JD  are  equal  at  the  corner.     This  approach  was  used  by  Lachat  and  Watson 


[27]  in  the  early  treatment  of  problems  in  3-d  elastostatics.  They  found  that  the 
resulting  errors  were  confined  mostly  to  the  region  of  the  domain  near  the  corner. 
Brebbia  and  Dominguez  [28]  also  investigated  this  method  with  regard  to  2-d 
potential  problems  but  found  that,  for  their  cases,  the  errors  in  the  gradients 
near  the  corners  were  unacceptable.  They  then  attempted  to  use  two  nodal 
points  in  close  proximity,  one  on  either  side  of  the  corner  (binodes),  concluding 
that  the  results  using  this  approach  were  generally  better  but  were  still  subject 
to  large  errors.  Riceurdella  [29]  used  essentially  the  saime  "binode"  approach  in 
his  treatment  of  problems  in  elasticity,  but  he  considered  only  step  changes  in 
tractions  (Case  3),  or  traction/displacement  combinations  (Cases  1  and  2), 
avoiding  Case  9  cdtogether. 


TABLE  1-2.  The  unknowns  and  corresponding  number  of 
equations  available  in  the  standard  RVBEMs  for 
the  nine  possible  corner  node  boundary  condition 
combinations  listed  in  Table  1-1. 


Number  of 
Case  Unknowns  Eqns.  Available 


1 

1 

2 

1 

3 

^ 

1 

4 

#.   and   i^). 

2 

5 

#.and(|^),          r 

-w^^- 

6 

^dr)^^    ^dr)^^    ^"^    '^ 

3 

7 

(qJd,   and  (f> 

2 

8 

ign)uy   and   <f) 

2 

9 

(^^^  -'  (i)^ 

1 

Alaxcon  et  al.  [30]  proposed  a  different  approach  based  on  the  fact  that,  in 

principle,  knowledge  of  the  potential  distribution   along  the  boundary  allows 

calculation  of  the  gradient  there.    By  cissuming  a  Unear  variation  of  -^  and  -^ 

on  OS 

over  the  segments  connecting  the  corner  node  with  its  adjacent  nodes,  each  of 
the  normal  gradients  is  expressed  in  terms  of  a  single  "new"  unknown,  V^. 
After  solution,  the  values  of  C^Xj  and  (o^)^)  are  calculated  from  V  <f).  The 
drawback  with  this  method  is  that  as  the  angle  turned  by  the  tangent  to  the 
boundary  decreases,  the  approximation  to  V  (^  becomes  less  and  less  accurate. 
Paris  et  al.  [31]  later  proposed  calculating  the  values  of  (o^h  and  (-^Jd  directly 
from  the  values  of  the  potential  at  the  corner  and  corner-adjacent  nodes.  This  is 
done  before  the  BEM  solution  process  is  even  begun.  The  equations  associated 
with  the  corner  nodes  are  thus  completely  eliminated  from  the  BEM  system  of 
equations;  however,  such  direct  approximation  of  the  gradients  is  subject  to 
truncation  errors  and,  as  before,  the  approximations  to  the  normal  gradients 
become  unreliable  as  the  corner  becomes  less  sharp.  Walker  and  Fenner  [32] 
carried  the  "gradient  approximation"  line  of  thought  in  a  slightly  different 
direction  by  combining  approximate  expressions  for  (-^\j  and  (-^jo  into  a  single 
equation  relating  the  two.  This  equation  is  then  used  as  the  "extra"  equation 
needed  to  properly  constrain  Case  9.  The  drawback  associated  with  this  method 
is  that  the  equation  relating  the  normal  gradients  is  sometimes  inaccurate,  with 
special  treatment  being  necessary  when  the  tangent  to  the  boundary  turns 
through  an  angle  neax  ninety  degrees  or  when  the  tangential  gradients  are  nearly 
zero  and  the  angle  turned  by  the  tangent  is  less  than  thirty  degrees. 

Grilli  et  al.  [33]  have  reported  good  results  using  a  similar  "gradient 
approximation"  approach  in  their  investigation  of  non-linear  waves,  as  have 
Bruch  and  Lejeune  [34]  in  their  study  of  2-  and  3-d  potential  problems.     The 


8 

latter's  treatment  is  somewhat  problem  dependent,  however,  since  they  state 
that  "the  additional  equation  (for  Case  9)  will  vary  from  problem  to  problem." 

Yet  another  approach  was  investigated  by  Gray  [35],  who  used  the  so- 
called  hypersingular  equation  (obtained  by  differentiating  the  traditional 
boundary  integral  equation)  as  the  "extra"  equation  needed  in  Case  9. 
Reasonable  results  were  obtained  for  practical  problems  involving  the  simulation 
of  a  3-d  electrochemical  plating  process.  The  drawback  associated  with  this 
method  is  the  careful  treatment  required  for  proper  evaluation  of  the  integr«ds  in 
the  hypersingular  equation. 

A  completely  different  line  of  reasoning  involves  the  use  of  discontinuous 
(or  non-conforming)  boundary  elements  at  a  corner.    A  discontinuous  element  is 
one  in  which  the  collocation  points  aie  removed  from  the  element  perimeter  and 
do  not  coincide  with  the  geometric  nodes  used  to  prescribe  the  domain  geometry. 
Using  this  approach,  a  collocation  point  is  placed  on  either  side  of  a  comer  with 
the  customary  single  geometric  node  still  located  at  the  corner  (see  Fig.  1-1). 
The  corner  geometry  is  thus  retained,  while  the  use  of  two  collocation  points 
permits  the  writing  of  two  boundary  element  nodal  equations.    This  allows  the 
normal   gradient   of  the   potential   to   be   double-valued   at    the   corner.      One 
drawback  associated  with  this  approach  is  that  since  no  collocation  node  actually 
exists  at  the  corner,  the  6  and  ^  values  must  be  interpolated  from  values  at  the 
interior  collocation  points.     This  introduces  some  inaccuracies  at  the  corners. 
Also,    the    number    of    simultaneous    equations    is    increased    by    the    use    of 
discontinuous  elements  since  two  nodes  are  used  near  each  corner  as  opposed  to 
one  node  for  continuous  elements.    Patterson  and  Sheikh  [36]  successfully  applied 
this  method  to  2-d  harmonic  and  3-d  potential  problems,  and  Danson  et  al.  [37] 
have   made   use   of  discontinuous   elements   in   their   general   purpose   BEASY 
boundary  element  program  with  good  results.     However,  Manolis  and  Banerjee 


^ — --p^    separate 

y^    /^    collocation  points 


ordinary 
node 


corner  geometric 
node 


FIGURE  1-1.   Using  "discontinuous"  elements  at  a  boundary  corner. 
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[38]  compared  the  use  of  continuous  (conforming)  and  discontinuous  (non- 
conforming) boundary  elements  and  found  that,  in  many  cases,  the  use  of 
discontinuous  elements  lead  to  large  errors  in  the  computed  solution.  They 
attributed  these  errors  to  the  fact  that  discontinuous  element  collocation  nodes 
often  do  not  lie  on  the  actual  domain  boundary,  concluding  that,  "In  general, 
conforming  (continuous)  elements  yield  more  accurate  results  than  non- 
conforming (discontinuous)  elements."  This  conclusion  was  subsequently 
challenged  by  Brebbia  and  Niku  [39]  who  suggested  that  the  errors  attributed  to 
the  use  of  discontinuous  elements  were  actually  the  result  of  poor  numerical 
methods — a  chcirge  subsequently  denied  by  Manolis  and  Banerjee.  It  appears 
that  the  use  of  discontinuous  elements  is  still  under  debate  in  BEM  circles; 
nevertheless,  the  BEASY  code  continues  to  use  them  with  success. 

Recently,  Detournay  [40]  has  investigated  the  handling  of  comer 
singularities  using  a  complex  variable  integral  equation  method  similar  to  the 
CVBEM.  In  that  paper,  special  corner  elements  characterized  by  two 
intersecting  straight  line  segments  were  used  to  handle  problems  where  the 
gradient  of  the  potential  had  a  constant  direction  along  each  of  the  segments. 
By  using  these  elements,  it  is  intended  that  the  special  behavior  of  the  potential 
function  in  the  vicinity  of  a  corner  can  be  more  accurately  modeled.  A  drawback 
with  this  approach  is  that,  in  many  cases,  such  an  element  yields  a  non-linear 
nodal  equation. 

By  formulating  2-d  potential  problems  with  complex  variables  instead  of 
real  variables,  the  CVBEM  avoids  the  problem  of  gradient  discontinuity  at 
boundary  corners.  Instead  of  solving  for  the  double-valued  gradients  at  a  corner 
node,  one  now  solves  for  a  single-valued  stream  function.  Still,  Neumann  and 
Robin  conditions  can  lead  to  special  circumstances  at  corners,  cind  a  methodology 
for  the  treatment  of  corners  in  the  CVBEM  is  needed. 
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1.2  Higher- Order  Elements 

The  earliest  BEM  formulations  approximated  the  domain  boundaxy  as  a 
series  of  linear  segments  over  which  the  value  of  the  variables  of  interest 
(temperature,  displacement,  etc.)  were  constant  [41,42].  These  segments  axe 
known  as  constant  elements  (see  Fig  l-2a).  The  collocation  points  where  the 
unknown  values  are  considered  in  the  BEM  solution  are  known  as  nodes  or  nodal 
points,  and,  for  constant  elements,  they  are  located  right  at  the  center  of  each 
element.  Constant  elements  provided  reasonable  solutions  in  many  cases,  but 
when  greater  accuracy  was  required,  investigators  used  more  refined  boundary 
approximations  such  as  linear,  quadratic,  or  cubic  elements  [1,2,28,29].  Linear 
elements  have  nodes  at  both  ends  (see  Fig  l-2b),  and  the  variables  of  interest  are 
assumed  to  vary  linearly  over  each  element.  Quadratic  and  cubic  elements  are 
referred  to  as  higher-order  elements,  and  they  allow  the  curved  structure  of 
boundaries  to  be  retained.  As  their  names  imply,  quadratic  elements  pass  a 
second-degree  polynomial  through  three  successive  nodal  points,  while  cubic 
elements  involve  third-degree  polynomials  and  four  nodal  points  (see  Figs.  l-2c 
and  l-2d).  .  . 

Hromadka  [7]  has  Iciid  the  groundwork  for  the  use  of  higher-order  elements 
in  the  CVBEM  via  generalized  proofs,  but  linear  elements  remain  the  highest- 
order  elements  that  have  been  reported  to  date  in  the  CVBEM  literature. 

1.3  Numerical  Grid  Generation 
It  was  mentioned  earlier  that  finite-difference  methods  (FDMs)  are  widely 
used  in  the  solution  of  PDEs,  especially  in  the  field  of  computational  fluid 
dynamics  (CFD)  [43,44].  All  FDMs  require  that  the  continuous  domain  of 
interest  be  replaced  by  a  discrete  domain  composed  of  a  collection  of  points 
distributed  throughout  the  original  continuous  domain.    These  points  are  known 
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FIGURE  1-2.    Different  types  of  boundary  elements. 
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as  grid  points,  and  the  collection  of  grid  points  is  known  as  a  grid  system. 
Proper  placement  of  these  grid  points  is  essential  if  accurate  solutions  axe  to  be 
obtained.  Often,  nonuniform  distributions  are  necessary  in  order  to 
acconmiodate  irregular  geometries  or  complex  flow  conditions.  Unsteady 
problems  may  even  require  grid  points  which  move  over  time.  Unfortunately,  it 
is  extremely  difficult  and  impractical  to  derive  finite-difference  equations  (FDEs) 
with  non-uniformly  distributed  and  moving  grid  points.  For  this  reason,  the  grid 
points  in  the  "physical"  spatial  domain  are  mapped  onto  a  "transformed" 
domain  where  they  are  stationary  and  uniformly  distributed  (see  Fig.  1-3).  In 
practice,  it  is  easier  to  perform  the  mapping  in  the  opposite  direction,  and  this 
process  is  known  as  grid  generation,  i.e.,  the  generation  of  grid  points  in  the 
"physical"  spatial  domain  [45]. 

Methods  for  generating  grid  systems  are  traditionally  divided  into  two 
major  classes — differential  equation  methods  (including  conformal  mappings)  and 
algebraic  methods.  Differential  equation  methods  generate  grids  by  solving  one 
or  more  PDEs  that  describe  the  transformation  between  the  "physical"  spatial 
domain  and  the  "transformed"  domain.  This  usually  requires  a  significant 
computational  effort,  but  it  produces  grid  systems  whose  grid  lines  are  smooth 
and  non-overlapping.  On  the  other  hand,  algebraic  methods  generate  grid 
systems  by  interpolating  between  the  boundaries  of  the  "physical"  spatial 
domain.  Since  the  solution  of  PDEs  is  not  involved,  the  algebraic  methods  tend 
to  be  much  more  computationally  efficient  than  the  differential  equation 
methods.  However,  grid  systems  generated  by  algebraic  methods  may  have  grid 
lines  which  are  non-smooth  and  overlapping — characteristics  which  must  be 
corrected  if  meaningful  finite-difference  solutions  are  to  be  obtained. 

Thompson  et  al.  [46]  have  written  a  book  which  thoroughly  describes  the 
techniques  used  in  numerical  grid  generation.     Also,  two  comprehensive  survey 
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papers  on  grid  generation  are  listed  as  References  47  and  48  for  the  interested 
reader.  Here,  only  those  grid  generation  methods  which  axe  in  some  way- 
connected  to  the  CVBEM  are  reviewed. 

Since  the  CVBEM  utilizes  complex  variables,  it  shares  common  ground 
with  the  conformed  mapping  techniques.  Generally,  these  techniques  have  been 
restricted  to  particular  geometries,  but  some  of  the  more  flexible  approaches 
include  those  of  Anderson  [49],  Davis  [50],  Halsey  [51],  and  Harrington  [52].  The 
CVBEM  is  set  apart  from  even  these  more  advanced  conformaJ  mapping 
methods  in  terms  of  its  ability  to  handle  irregular  geometries  due  to  its  boundary 
element  nature.  For  this  reason,  references  which  involve  the  application  of 
BEMs  to  grid  generation  will  now  be  reviewed. 

The  integral  equation  method  of  Symm  [42],  mentioned  earlier,  has  been 
used  for  mapping  a  simply-connected  2-d  spatial  domain  onto  a  unit  disk  and  for 
mapping  doubly- connected  region  onto  annular  regions.  Hayes  et  al.  [53] 
subsequently  improved  this  approach  and  were  able  to  obtain  more  accurate 
results.  In  each  of  these  applications,  the  primary  concern  was  to  obtain  the 
mapping  rather  than  to  use  the  mapping  for  generating  grid  systems.  Indeed, 
the  only  application  of  a  BEM  directly  to  grid  generation  appears  to  be  that 
described  by  Adamczyk  [54],  who  used  an  electrostatic  analog  to  generate  grids 
between  bodies  in  a  cascade.  Force  and  potential  lines  calculated  using  a  panel 
method  solution  of  an  integral  equation  were  used  to  define  the  resulting  grid 
system.  Also,  a  doubly  infinite  line  of  alternating  charges  separated  by  a  finite 
distance  was  used  as  the  fundamental  solution. 

It  appears  that  the  CVBEM  has  never  been  applied  to  numerical  grid 
generation;  however,  its  accuracy  and  efficiency  make  it  an  excellent  candidate 
for  bridging  the  gap  between  existing  algebraic  and  differential  equation 
methods.     Since  it  involves  the  solution  of  the  Laplace  equation,  it  would  be 
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properly  classified  as  a  differential  equation  method  (although  the  term  integrcd 
equation  method  would  be  more  appropriate),  but  it  is  expected  to  be  much 
faster  than  the  existing  elliptic  solvers. 

1.4  Objectives 

Based   on   the   background   information   presented   in   the  previous   three 
sections,  the  objectives  of  this  dissertation  are  as  follows: 

1)  Develop  a  formulation  of  the  CVBEM  using  quadratic  elements. 
Previously,  lineax  elements  were  the  highest-order  elements  used. 

2)  Devise  a  methodology  for  the  proper  treatment  of  boundary  corners  in 
the  CVBEM.  This  topic  has  been  treated  extensively  in  the  RVBEMs 
but  has  not  been  addressed  in  the  CVBEM. 

3)  Derive  a  variation  of  the  CVBEM  which  can  be  used  for  numerical  grid 
generation.  It  is  noted  that  the  CVBEM  has  never  been  used  for  this 
purpose. 

1.5  Overview 

In  the  next  chapter,  a  description  of  the  linear  element  CVBEM  as 
formulated  by  Hromadka  is  presented  together  with  some  observations  and 
suggestions  regarding  some  of  the  method's  more  subtle  points.  Chapter  III  then 
extends  the  CVBEM  via  the  use  of  quadratic  elements.  Chapter  IV  goes  on  to 
discuss  the  proper  treatment  of  corners  in  both  the  linear-  and  quadratic-element 
CVBEMs,  while  Chapter  V  presents  a  novel  application  of  the  CVBEM  in  the 
area  of  numerical  grid  generation.  Examples  and  results  which  demonstrate  the 
new  techniques  of  the  previous  three  chapters  are  contained  in  Chapter  VI. 
Finally,  in  Chapter  VII,  conclusions  are  stated  and  recommendations  for  further 
work  are  given. 
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To  present  the  material  in  a  concise  manner,  all  derivations  are  relegated 
to  the  Appendix. 


CHAPTER  II 
DESCRIPTION  OF  THE  LINEAR-ELEMENT  CVBEM 


When  solving  2-d  potential  problems,  it  is  often  convenient  to  make  use  of 
the  theory  of  complex  variables.  Using  this  approach,  one  can  construct  a 
complex  potential,  uj{z)  =  <f>{x^y)  +  i0(x,y),  where  (j)  is  the  potential  function 
and  ^  is  the  stream  function.  This  complex  potential  is  analytic  so  that  (f)  and  ip 
satisfy  the  Cauchy-Riemajin  conditions, 


where  s  represents  the  streamline  coordinate  and  n  represents  the  coordinate 
normal  to  it;  furthermore,  one  can  apply  the  well-known  Cauchy  integral  formula 
which  states  that  <    "       '  '■■--,. 

r 

This  expression  relates  the  value  of  complex  potential  u  at  point  Zq  located 
inside  the  k-connected  Jordan  domain,  fi,  to  a  contour  integral  (containing  w) 
along  the  boundary,  F.  The  direction  of  travel  for  the  contour  integral  is  such 
that  the  interior  of  the  domain  is  always  to  the  left. 

The  linear-element  formulation  of  the  CVBEM  transforms  the  Cauchy 
integral  formula  into  a  BEM  by  using  two  major  approximations.  First, 
boimdary  F  is  discretized  into  N  finite- length  segments  (elements),  Tj,  the  end 
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points  of  which  axe  referred  to  as  nodal  points  or  nodes.    Domain  boundary  F  is 
then  taken  as  the  union  of  these  elements  as  shown  in  Fig.  2-1,  i.e., 


r=u  r,  (2.3) 


Second,  the  function  ui{z)  is  approximated  by  a  lineax  global  trial  function, 
Gi(z),  given  by 


Gi(^)   =    E  N,(^)  u;,-  (2.4) 

>  =  1 


where  uij  is  the  value  of  u  at  nodaJ  point  Zj,  and  ^j{z)  is  a  continuous  basis 
fimction  weighting  the  effect  of  ufj  over  elements  Tj  _  j  and  Tj. 

It  is  necessary  to  explain  the  notation  that  wiU  be  used  throughout  the 
remainder  of  this  dissertation.  The  subscript  exact  will  identify  a  quantity  whose 
value  is  known  exactly.  The  bar,  ~ ,  will  refer  to  a  quantity  whose  value  is 
specified,  such  as  in  the  boundary  conditions.  Finally,  the  hat,  "",  will  represent 
a  quantity  whose  value  is  treated  as  an  unknown  in  the  solution  by  the  CVBEM. 

Returning  now  to  the  description  of  the  linear  element  CVBEM,  basis 
function  ^j{z)  is  taken  as  a  first  degree  polynomial  of  the  following  form: 

'  {z-zj_,)  /  {zj-zj_^)  z  e  r^_i 

Ni(^)     =      I      (^i  +  l-^)/(2i  +  l-2j)  Z     e     Tj 

0  z  ^r,._,  u  r, 


By  substituting   Gi{()  for  uj(Q  in  the  right  hand  side  of  the  Cauchy 
integral  formula  (Eqn.  (2.2)),  the  first-order  CVBEM  approximation  of  u  can  be 
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e(j+i,j-i;j 


FIGURE  2-1.      Boundajy  discretization  and  angle  definitions  in  the  lineax- 
element  CVBEM. 
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Gi(C). 


^(^o)   =  i  /  r^C 


zq  e  ii     ^0  ^  r 


(2.5) 


After  much  simplification  (see  [7]),  the  contour  integral  can  be  eliminated,  and 
the  above  equation  reduces  to 

.        N  ft. 

^^^°^   "  2ii  5Z  K  +  i(^o  -  ^j)  -  ^ji^o  -  Zj^x)\-r ^-— T  (2.6a) 


where     h^  =  In 


(^i  +  i-^o) 


(^j  -  2o) 


=  In 


(^j  +  i-^o) 


(^i  -  ^o) 


+  t-5(j  +  l,  i;0)       (2.6b) 


This  formula  forms  the  basis  for  the  linear-element  CVBEM.     (See  Fig.  2-1  for 
the  definition  of  ^(7 -I- 1,  j;  0).) 

Eqn.  (2.6),  being  expressed  in  complex  variables,  actucilly  embodies  two 
equations — one  for  the  real  part  and  one  for  the  imaginary  part.  If  the  values  of 
^  and  t/)  (and  thus  w)  are  known  at  each  boundary  node,  Eqn.  (2.6)  can  be  used 
to  calculate  <^  eind  V'  at  any  interior  point,  Zq.  In  most  potential  problems, 
however,  boundary  conditions  specify  either  <^  ,  ^,  or  neither  of  them  explicitly. 
To  solve  for  the  unknown  values  of  (^  eind  V',  it  is  necessary  to  derive  an  extended 
version  of  Eqn.  (2.6)  by  moving  Zq  to  the  position  of  Zj  on  the  boundary.  In  this 
effort,  one  cannot  simply  plug  in  Zj  for  Zq  in  Eqn.  (2.6)  because  {zj  —  Zq)  appears 
in  the  denominator  of  the  natural  log  term.  Instead,  one  takes  the  limit  of  Eqn. 
(2.5)  as  Zq  approaches  Zj.  Hromadka  [7]  has  performed  this  somewhat  involved 
task,  with  the  result 


u;,  =  TT^  Wi 


In 


^3  +  1        ^3 


'3-l~^3 
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+  eu  +  h  i-i;  i) 


where   h-  =  In 


N  '  1 

i^i  +  l-^j) 


2.  +  l-^i 


z-  —  z  ■ 
«        J 


=  In 


(^.-  -  ^j) 


+  ie{i  +  l,  i;  j) 


(2.7a) 


(2.7b) 


The  angles  6{j  +  1,  ji  —  1;  j)  and  ^(i  +  1,  i;  j)  axe  shown  in  Fig.  2-1. 

Eqn.  (2.7)  can  be  applied  at  any  boundaxy  node,  but  like  Eqn  (2.6),  Eqn. 
(2.7)  has  real  and  imaginary  parts.  Two  equations  can  thus  be  derived  for 
arbitrary  boundary  node  j  as 


1 


'f>j=  2^    V'.ln 


^i  +  i-^j 


2j  _  1  -  Zj 


+    <l>.0{j  +  l,   j-l-j) 


+  E  {<?^.+iC2  +  v-.+iCi  -  <^,C4  -  tz-A}} 


•.•  +  1#J 


and 


'^^=    -2^V^i^^ 


^3  +  l~  ^3 


^3-1-^3 


-  4,^e{j  +  l,j-l;j) 


'    ! 


+  E  {<^.  +  iCi  -  ^.  +  102  -  <A.C3  +  v-.c, 


where 


I  =  1 


Ci  =  [ix,-x,)C  -  {yj-vM 

C2  =  [(x,-x,)D  +  (y,-j/.)C] 

C3  =  [(xj-Xi^,)C  -  (y,-y.  +  i)D] 

C4  =  [(x,-x,  +  i)D  +  (j/,-j/i  +  i)C] 

C  =  [A(x,  +  i-x.)  +  B(y..,i-y,)]/F 


(2.8) 


(2.9) 

(2.10a) 
(2.10b) 
(2.10c) 
(2.10d) 
(2.10e) 
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D  =  [B(x.  ^  1  -  X.)  -  A(y.  + 1  -  y.)]  /  F  (2.10f) 


F  =  (x,  +  ,-x,)2  +  {y,^,-y.f  (2.10g) 

A  =  ln|^^^|  (2.10h) 

B  =  9{i  +  l,i;  j)  (2.10i) 


Henceforth,  Eqns.  (2.8)  and  (2.9)  will  be  referred  to  as  the  phi  nodal  and  psi 
nodal  equations,  respectively. 

2.1  Boundary  Conditions 
Having  generated  the  desired  equations  for  calculating  (f)  and  ^  at  the 
boundary  nodes,  a  question  arises  as  to  how  one  uses  Eqns.  (2.8)  and  (2.9)  to 
solve  for  the  unknown  boundary  values  of  (f>  and  0.  Before  this  question  can  be 
answered,  consideration  must  be  given  to  the  conditions  that  can  exist  at  the 
boundary  nodes.  As  it  turns  out,  each  of  the  three  "classic"  boundary  conditions 
of  the  potential  theory  yields  a  different  equation  as  detailed  below  [21]. 

Dirichlet  Condition:    The  potential  at  node  _;  is  known  and  specified,  i.e., 

?,    =    Kact.j  '  '      '       .       "     '^  ■•        (2.11) 

Notice  that  for  such  a  condition  0j  is  unknown.  In  keeping  with  the 
notation  defined  previously,  the  specified  <f)j  becomes  ^j,  and  the  unknown 
ipj  becomes  V'j- 

Neumann  Condition:  The  normal  gradient  of  the  potential  at  node  j  is 
known.    The  stream  function  there  is  then  related  to  the  normal  gradient 
I  of  the  potential  as 
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^^-'^i-  +H©,  +  0.-.}K-S-.l  (2.12) 

Here,  it  is  cissumed  that  (^j  varies  linearly  from  node  j  —  1  to  node  j  [21]. 
For  this  condition,  both  (^j  and  tpj  axe  unknown  and  aire  referred  to  as  <^  • 
and  t/jj. 


Robin    Condition:        The    normal    gradient    of    the    potential    and    the 
potential  itself  axe  related  by  -^ 
is  then  related  to  the  potential  as 


potential  itself  are  related  by  ^  =   H(<^  -  <f>g^).    The  stream  function  there 


^,-  =  V-,- _  1  +  ^  {  (H(^^  - ^  +  (H(<^^  -  <j>)l _ ,}  \z^  -z,_,\  (2.13) 

Here,  it  is  assumed  that  H<^  and  H^qq  vary  linearly  from  node  j  —  \  to  node 
j  [21].  As  in  the  Neumann  condition,  both  4>j  and  V'>  are  unknown  and 
again  become  ^ J  and  V* J. 

The  use  of  complex  variables  in  the  CVBEM  also  allows  for  the  occurrence 

of  the  stream  function  condition  as  follows: 

"i 

Stream  Function  Condition:  The  stream  function  at  node  j  is  known  and 
specified,  i.e.,  ^ 

^J    =    V'.xact,;  (2.14) 

For  this  condition,  4>j  is  unknown.    The  specified  V'j  is  renamed  V'j,  and  the 
unknown  (j)j  becomes  (j)j. 
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2.2  Solution  Methods 
Eqns.  (2.8)  and  (2.9)  will  be  used  to  estimate  the  unknown  values  of  (f>  and 
%p  at  the  boundary  nodes,  but  prior  to  the  presentation  of  the  methodology  for 
doing  so,  a  close  examination  of  these  equations  is  in  order.  It  is  clear  from  their 
format  that  (j>j  and  V'j  on  the  left  can  be  calculated  by  using  the  known  values  of 
<j)  and  V'  S't  all  boundary  nodes  whose  indices  appear  on  the  right.  As  discussed 
previously,  some  of  the  vcdues  of  (f)  and  0  at  the  boundary  nodes  are  not  specified 
by  the  boundary  conditions.  The  methods  for  estimating  the  unspecified  values 
of  <f>  and  xl>  (designated  (f)  and  t/))  thus  hinge  on  how  these  quantities  are  related 
to  the  specified  qucintities  {(f)  and  %l>)  and  on  which  of  the  equations  ((2.8)  or 
(2.9))  is  used  in  the  construction  of  the  matrix  equation  for  solution  of  the 
problem.  At  nodal  point  j  where  a  Dirichlet  or  stream  function  condition  is 
imposed,  there  are  three  methods  of  solving  for  the  estimated  nodal  values  of  <i)j 
or  V'j. 

Explicit  Method:  For  a  Dirichlet  condition  imposed  at  node  j,  <j)j  is 
known  (as  (l>j)  but  ^^  is  not.  One  thus  sets  (f>j  =  (l>j  =  (f>j,  and  rpj  =  xl)j.  The 
first  setting  is  governed  by  the  fact  that  a  Dirichlet  condition  is  specified; 
no  estimation  is  thus  needed  for  (f)y  The  second  setting  is  made  in  order  to 
estimate  the  unknown  value  of  V'j-  Without  such  a  setting,  there  will  be 
two  unknowns  (V'j  and  0^)  at  the  same  nodal  point,  a  situation  which  does 
not  allow  for  a  unique  solution.  Since  the  field  theory  predicts  that  the 
estimated  V'j,  if  error  free,  will  be  exactly  equal  to  that  imposed,  equating 
xj^j  and  %l)j  is  certainly  reasonable. 

In  the  explicit  method,  an  effort  is  made  to  keep  all  of  the  unknowns 
on  the  right  side  of  the  equation.  It  is  therefore  impossible  to  use  Eqn. 
(2.9)   since  ^j  on   the  left   side  is  unknown  for  the  Dirichlet   condition 
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specified.    Thus,  Eqn.  (2.8)  is  used  as 


1    )j. 


Zi,y    -  Z, 


+   E  {<^.  +  iC2  +  V.  +  iCi  -  <i>,C,  -  0.C3}}  (2.15) 


t  =  1 

«.«  +  !  94  J 


This  equation  is  referred  to  as  the  explicit  phi  nodal  equation.  For  the 
explicit  method,  it  is  sufficient  to  use  Eqn  (2.15)  to  solve  for  ^j.  Eqn.  (2.9) 
is  dropped. 

For  a  stream  function  condition  at  node  j,  xpj  is  known  (as  ipj)  but  ^ 
is  not.    One  thus  sets  0^  =  xl)-  =  ipj,  and  (f>j  =  Zj.    This  time,  for  the  expHcit 
method  of  solution,  Eqn  (2.9)  is  used,  and  the  nodal  equation  is 


^,=  -^{?,la 


^i  + 1      ^i 


^j-i-^j 


-  V,-^(i  +  l,  j-l;j) 


.    \       +  E  {<^.  +  iCi  -  0.  +  1C2  -  <?!>.C3  +  0,04}!  (2.16) 

'  '         '•'+^^J  -  .  ..     .       .     .  .,_. 

This  equation  is  referred  to  as  the  explicit  psi  nodal  equation.    Eqn.  (2.8)  is 
dropped  for  this  node  in  the  solution. 

For  the  explicit  method,  solution  is  effected  by  applying  Eqn.  (2.15) 
at  each  Dirichlet  node  and  Eqn.  (2.16)  at  each  stream  function  node  to 
form  a  system  of  simultaneous  equations  which  is  solved  for  the  unknown 
values  of  (^  and  xj). 

Implicit  Method:  For  a  Dirichlet  condition  imposed  at  node  j,  one  sets 
(i>j  =  (f>j,  and  tj^j  =  xj^j.  However,  unHke  the  explicit  method,  the  impHcit 
method  involves  a  nodal  equation  where  unknowns  appears  on  both  sides. 
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Eqn.  (2.9)  is  thus  used  as 


^>=  "ij-^^^^ 


^  j + 1   ^j 


^j-i-^j 


N 


-  4,.0ij  +  l,j-l;j) 


+   E  {<^.  +  iCi  -  V'.  +  iC^  -  .^A  +  ^Ajj        (2.17) 

This  equation  is  referred  to  as  the  implicit  psi  nodaJ  equation.    Eqn.  (2.8) 
is  dropped. 

Similaxly,  for  a  stream  function  condition  specified  at  node  j,  one  sets 

tpj  =  i/^j  and  (^j  =  (f)j.   Eqn.  (2.8)  is  now  used  as 


1  h. 


^i  + 1      ^3 


^>  ^  s^FV^i^^  C\-z]  +  <?^.%'  +  i.  i-i;j) 


+  E  {<;^.-+iC2  +  tA.+iCi  -  <i>,c,  -  r^ 


I  =  1 


A}} 


(2.18) 


This  equation  is  referred  to  as  the  implicit  phi  nodal  equation.    Eqn.  (2.9) 
is  dropped.  ,     ..„  ^ 

In  a  maimer  similar  to  the  explicit  method,  solution  by  the  implicit 
method  is  effected  by  applying  Eqn.  (2.17)  at  each  Dirichlet  node  and  Eqn. 
(2.18)  at  each  stream  function  node  to  form  a  system  of  simultaneous 
equations  which  is  solved  for  the  unknown  values  of  (j>  and  tj). 


Hybrid  Method:  For  a  Dirichlet  condition  imposed  at  node  j,  one  sets 
(f>j  =  (j)j  on  the  left  side  of  Eqn.  (2.8)  aoid  (f>j  =  (j)j  and  V';  =  '^j  on  the  right 
sides  of  Eqns.  (2.8)  and  (2.9)  to  obtain 
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_    1 


^  j  + 1     ^] 


•^i  =  2^  V^^- ^^  ^TTT^    +  .^,-^(;  +  l,  ;-l;i) 


+   L  {<?^.  +  iC2  +  V'.  +  A  -  <i>,C,  -  0 


i  =  1 
«.•  +  !  #i 


A}} 


(2.19) 


0,=    - 


^{^.1- 


^j  +  i-2^j 


^j-i-^j 


-    0.0(j  +  l,   j-l;j) 


+  E  {^.+iCi  -  V'.+iCj  -  ^A  +  v-a}} 


>, '  + 1  #  i 


(2.20) 


These  equations  axe  referred  to  as  the  Dirichlet  hybrid  phi  nodal  and  psi 
nodal  equations,  respectively,  since  Eqn.  (2.19)  is  explicit  while  Eqn.  (2.20) 
is  implicit.    Notice  that  both  (t>j  and  ipj  are  unknown. 

For  a  stream  function  condition  imposed  at  node  j,  one  sets  0^  =  V', 
on  the  left  side  of  Eqn.  (2.9)  and  V*,  =  V'j  and  <^>  =  <l>i  on  the  right  sides  of 
Eqns.  (2.8)  and  (2.9)  to  yield 


.*        ■'  .'• 


^^■  =  i^.- 


In 


^3  + 1      ^i 


+  ^^(;  +  l,  i-l;i)        •    '■ 


+  E{4>i  +  iC,  +  0..,iCi  -  <^,C4  -  0A 


I  =  1 


(2.21) 


0,  =    -ifeln 


•^  j  + 1     ^j 


'j-i-^j 


-  0,-^(i  +  l,  i-l;i) 


+   E  {?^.  +  iCi  -  0i  +  iC2  -  <?i,C3  +  VA})  (2.22) 


•  =  1 
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These  equations  are  referred  to  as  the  stream  function  hybrid  phi  nodal 
and  psi  nodal  equations,  respectively,  since  Eqn.  (2.21)  is  implicit  while 
Eqn.  (2.22)  is  explicit.    Again,  both  4>j  aaid  0^  are  unknown. 

In  the  hybrid  method,  solution  is  effected  by  applying  Eqns.  (2.19) 
and  (2.20)  at  each  Dirichlet  node  and  Eqns.  (2.21)  and  (2.22)  at  each 
stream  function  node.  The  number  of  equations  solved  simultaneously  for 
unknowns  4>  and  tp  is  thus  doubled  in  the  hybrid  method  as  compared  to 
the  explicit  and  implicit  methods. 

Table  2-1  gives  a  summary  of  the  solution  methods  described  above  based 
on  a  node  imposed  with  a  Dirichlet  condition.  A  close  exeunination  of  this  table 
reveals  the  strengths  and  weaknesses  associated  with  each  method. 

Attention  is  first  directed  to  the  column  headed  "Quantities  Equated"  in 
the  table.  It  appears  that  the  only  seemingly  minor  difference  between  the 
explicit  method  and  the  imphcit  method  is  the  relation  between  (f>j  and  <^j,  i.e., 
the  two  are  set  equal  in  the  explicit  method  but  not  in  the  implicit  method.  In 
fact,  Eqn.  (2.8)  for  calculating  (f)j  is  dropped  in  the  implicit  method  when  solving 
for  0j  (see  the  colimin  headed  "Equation  Dropped").  However,  this  equation  can 
still  be  used  to  advantage  for  estimating  the  overall  accuracy  of  the  solution.  In 
this  effort,  (l)j  is  calculated  using  Eqn.  (2.8),  and  the  estimated  (^^  is  compared 
with  the  specified  (f)j  for  error.  A  "small"  error  gives  an  indication  of  high 
accuracy  in  the  overall  solution. 

The   above   advantage   does   not   apply   to  the  explicit   method.      For  a 

Dirichlet  condition,  the  explicit  method  requires  that  Eqn.  (2.9)  be  dropped.    If 

after  solution,  Eqn.  (2.9)  is  subsequently  used  to  estimate  the  unknown  0,,  no 

specified  V'j  exists  to  compare  it  to.     Worse  still,  errors  in  V'  at  the  boundary 

nodes   lead   to  errors   in   the   calculated   values   of  -t^,   since  -r^  is  found  bv 

on  an  ^ 
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TABLE    2-1.      Summaxy    of    solution    methods    based    on    a    node    imposed 
with  a  Dirichlet  condition. 


Method    Known  Unknown 


Quantities 
Equated 


Equation 
Used 


Equation 
Dropped 


Explicit 

^i 

^i 

<h  =  <t>j  =  <f>eTact,j 

(2.8) 

(2.9) 

Implicit 

h 

^i 

<t>j  =   <t>tXCiCt,j 

(2.9) 

(2.8) 

Hybrid 

^i 

V'i 

^j  =  ^j  =  <l>exact,i 

on  the  left  side 
of  Eqn.  (2.8); 

(2.8)  and  (2.9) 

none 

j   \ 

■ . 

-     <l>j  =  (t>j  and  rpj  =  V'j 
on  the  right  side  of 
•  Eqns.  (2.8)  and  (2.9) 
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numerically  differentiating  ip  as  will  be  discussed  later.  Experience  has  shown 
that  the  implicit  method  generally  yields  more  accurate  values  of  xp  than  the 
explicit  method. 

The  hybrid  method  is  essentially  the  combination  of  the  explicit  and 
implicit  methods.  As  such,  one  might  expect  that  it  is  the  m.ost  accurate  of  the 
three  methods  developed,  and  this  often  (but  not  always)  turns  out  to  be  the 
Ccise.    Table  2-1  shows  that,  in  the  hybrid  method 

(f>j  =  4)j  on  the  left  side  of  Eqn.  (2.8); 

(pj  =  ^j  on  the  right  sides  of  Eqns.  (2.8)  and  (2.9); 

tpj  =  V'j  on  the  right  sides  of  Eqns.  (2.8)  and  (2.9). 

The  fact  that  4>j  is  equated  to  (pj  on  the  left  side  of  the  Dirichlet  hybrid  phi 
nodal  equation  (Eqn.  (2.19))  renders  this  equation  very  similar  to  Eqn.  (2.15)  in 
the  explicit  method.  Likewise,  the  Dirichlet  hybrid  psi  nodal  equation  (Eqn. 
(2.20))  is  much  like  Eqn.  (2.17)  in  the  implicit  method.  The  differences  He  in 
the  fact  that  unknowns  (j)j  and  xpj  are  both  present  in  the  hybrid  method,  with 
Eqns.  (2.19)  and  (2.20)  both  being  applied  at  Dirichlet  node  j.  Since  ^j  is  also 
specified,  it  may  be  compared  after  solution  with  the  estimated  <f>j,  thus  offering 
the  potential  for  error  detection  afforded  by  the  implicit  method.  Unfortunately, 
as  just  noted,  the  hybrid  method  requires  the  application  of  two  equations  at 
each  Dirichlet  or  stream  function  node,  while  the  explicit  and  implicit  methods 
require  only  one.  This  results  in  a  doubling  of  the  number  of  simultaneous 
equations  associated  with  Dirichlet  and  stream  function  nodes. 

Be  mindful  that  all  of  the  above  comments  concerning  the  explicit, 
impUcit,  and  hybrid  methods  have  been  made  for  the  Dirichlet  case;  however, 
they  are  just  as  valid  for  the  stream  function  case  with  the  exchange  of  4>  with  ip. 
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For  nodes  where  Neumann  or  Robin  conditions  axe  specified,  Eqn.  (2.12)  or 
Eqn.  (2.13)  will  be  used  in  conjunction  with  the  implicit  phi  nodal  equation, 
Eqn.  (2.18),  in  the  solution  of  unknowns. 

By  applying  Eqns.  (2.12)  and  (2.18)  at  all  Neumann  nodes,  Eqns.  (2.13) 
and  (2.18)  at  all  Robin  nodes,  and  either  the  explicit,  implicit,  or  hybrid  method 
equations  at  all  Dirichlet  and  stream  function  nodes,  a  system  of  lineax  algebraic 
equations  is  obtained.    This  system  can  be  represented  in  a  matrix  form  as 


(II  ■« 


(2.23) 


Here,  C  is  the  square  matrix  of  coefficients  on  the  unknown  values  of  (f)  and  V" 
and  R  is  a  vector  of  known  constants.  A  detailed  description  of  how  C  and  R  axe 
formed  is  given  in  Appendix  A.  The  system  can  be  solved  by  any  direct  method 
such  as  Gaussian  elimination.  Once  the  unknown  values  of  (f>  and  V'  are  found  at 
each  boundary  node,  one  can  use  these  values,  together  with  the  specified  values, 
<i>  and  xj),  as  the  boundary  nodal  values  needed  by  Eqn.  (2.6)  for  calculating  Q  at 
any  interior  point.  ^      -■  ^     k  ~        .  ' 

The  matrix  formulation  of  Eqn.  (2.23)  provides  another  perspective  as  to 
the  preference  of  the  implicit  method  over  the  explicit  and  hybrid  methods.  In 
the  explicit  method,  the  nodal  equations  produce  a  coefficient  matrix  with  zeros 
on  the  main  diagonal — a  very  poor  situation,  numerically.  By  contrast,  the 
implicit  method  produces  a  matrix  with  the  largest  element  in  each  row  located 
on  the  main  diagoncd,  a  situation  conducive  to  obtaining  good  numerical  results. 
The  hybrid  method  suffers  from  the  same  numerical  pitfall  as  the  explicit 
method  (zeros  on  the  main  diagonal)  but  to  a  lesser  degree  since  implicit 
equations  are  also  used  to  form  the  hybrid  coefficient  matrix.  The  primaxy 
reason    for    choosing    the    implicit    method    over    the    hybrid    method    is    the 
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undesirable  increase  in  matrix  size  in  the  hybrid  method  due  to  the  need  for  two 
equations  at  each  Dirichlet  or  stream  function  node  instead  of  just  one. 

2.3  Calculation  of  the  Normal  Gradient  of  the  Potential 


After  ^  and  il>  have  been  evaluated  at  the  boundajy  nodes  and  interior 
s,  the  values  of  V"  may  be  used  to  deter 
is  made  of  the  Cauchy-Riemann  conditions, 


points,  the  values  of  V"  may  be  used  to  determine  ^  at  the  boundary  nodes.   Use 


d<i>  _  dip 


where  n  is  now  taicen  as  the  outward  normal  to  F  with  s  as  the  tangential 
coordinate  along  F,  assumed  positive  in  the  direction  of  the  contour  integral  in 
Eqn.  (2.2).  The  partial  differential  -^  (and  thus  ^)  can  be  approximated  by 
using  finite-difference  formulae.  For  example,  a  second-order  accurate,  three- 
point  backward  finite-difference  approximation  to  -J^  at  node  j  is  given  by 


ds  >'  (5_,  -  Sj  _  i)(5_,-  -  Sj  _  2)(5_,-  _  1  -  5_,-  _  2) 


j>i-j{Sj-Sj_^f 


{Sj  -  Sj  _i){sj  -  Sj  _  2){Sj  _  1  -  Sj  _  2) 


(2.24) 


The  values  of  the  tangential  axe  length,  s,  can  be  estimated  as  the  cumulative 
sum  of  the  lengths  of  the  linear  segments  connecting  each  successive  nodal  point, 
or  a  more  sophisticated  polynomial  or  spline  procedure  may  be  used.  One-sided 
formulae  should  be  used  at  corner  nodes  and  at  any  others  where  the  normal 
gradients  are  different  on  either  side  of  the  node.  Central-differencing  can  be 
used  at  all  remaining  nodes.    Equations  such  as  (2.24)  can  be  greatly  simplified  if 
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the  nodal  points  are  equally  spaced  along  the  boundary,  and  References  43  and 
44  provide  lists  of  some  common  fmite-difference  formulae  which  can  be  used 
when  this  is  the  Ccise. 

2.4  Additional  Comments  and  the  Psi  Reference  Node 
Two  more  points  are  worthy  of  note.  First,  nodal  points  which  are  located 
at  corners  can  be  considered  as  a  special  case  which  will  be  studied  in  great 
detail  in  Chapter  IV.  Second,  no  matter  what  boundary  conditions  exist,  one 
must  specify  a  reference  value  of  V"  at  some  nodal  point  along  the  boundary  in 
order  to  serve  as  the  constant  of  integration  in  Eqn.  (2.5).  Numerically  speaking, 
the  associated  matrices  become  singular  if  no  xj}  value  is  provided.  This  anchor 
node  will  be  referred  to  as  the  "psi  reference  node".  Numbering  it  as  node  ;,  one 
has 

r  5  ..  .  •  -  _ 

The  implicit  phi  nodal  equation  (Eqn.  (2.18))  can  then  be  used  to  determine  the 
unknown  (j)j  as  was  discussed  earlier  in  the  description  of  the  implicit  method. 

If,  however,  the  value  of  <t>exact  is  also  known  at  node  j  (i.e.,  a  Dirichlet 
condition),  one  can  specify  both  (f)  and  V'  at  this  node  as 

4>i     =    <Pexact,3 

and 

V'i    =    inexact,  j 

The  complex  potential  u  is  thus  fully  specified,  and  no  nodal  equation  needs  to 
be  applied.    This  node  will  be  called  a  "completely  specified  node". 
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In  some  circumstances,  the  position  of  the  psi  reference  node  caji  be 
important  to  the  application  of  the  CVBEM.  To  properly  discuss  this 
phenomenon,  a  detailed  examination  of  the  way  that  Neumann  and  Robin 
boundary  conditions  are  implemented  in  the  CVBEM  is  necessary.  Since  such 
an  examination  is  also  necesseiry  when  considering  treatment  of  comers,  further 
discussion  concerning  placement  of  the  psi  reference  node  will  be  delayed  until 
Chapter  IV  where  the  treatment  of  corners  is  covered. 

This  chapter  has  laid  the  foundation  for  the  development  of  the  CVBEM. 
A  linear-element  model  has  been  presented  here;  a  more  refined  model  follows  in 
the  next  chapter. 


CHAPTER  III 
FORMULATION  OF  THE  CVBEM  USING  QUADRATIC  ELEMENTS 


As  with  the  linear-element   CVBEM,  the  formulation  of  the  quadratic- 
element  CVBEM  begins  with  Cauchy's  integral  formula, 

r    . 

This  time,  however,  boundary  F  is  discretized  into  M  finite-length  curved 
elements,  F^.  These  elements  axe  formed  by  passing  a  quadratic  polynomial 
through  three  successive  nodaJ  points.  The  domain  boundary  is  taken  as  the 
union  of  these  elements  as  shown  in  Fig.  3-1,  i.e.. 


r=.U  r,       .      :  ^       ,    ^    k    .  !    .  (3.2) 

J  =  1 


Notice  that  the  element  _;'  contains  nodes  k,  k  +  1,  and  fc-i-2,  so  that  the  indices  j 

and  k  are  related  by  k  =  2j  -1.    Further,  A/  elements  are  formed  from  N  =  2M 

-,•■■'■"' 
nodal  points. 

The    function    u>{z)    is    now    approximated    by    a    quadratic    global    trial 

function,  G2(z),  given  by 


M 

G^{z)   =    '£Mk{z)uk  +  ^k  +  i{z)<^k  +  i  +  ^k  +  2{z)^k  +  2  (3.3) 

k  =  2j-l 
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8(k+2,k-2;k) 


ly 


-x 


FIGURE  3-1.      Boundary     discretization     and     angle     definitions    in     the 
quadratic-element  CVBEM. 
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Here,  Mfc(z),  Mjt  +  i(2),  and  M^^jl-^)  ^^^  continuous  basis  functions  representing 
the  effects  of  Uf.,  uji.^i,  and  uJk  +  2  o"^^r  elements  Tj-j  and  Tj.  These  basis 
functions  are  taken  as  second-degree  Lagrange  polynomials  of  the  form  \  \.  . 


M,iz)   = 


(^-^fc-2)(^-^fc-l) 
(^fc  +  l  ~  ^)  (^fc  +  2  ~  ^) 


•2^  e  r,_i 


z  e  r, 


z  ^r,_i  u  r,. 


(3.4) 


By  substituting  G2(C)  for  a?(C)  in  the  right  hand  side  of  the  Cauchy  integral 
formula  (Eqn.  (3.1)),  the  second-order  CVBEM  approximation  of  u  can  be 
expressed  as  .'  ■ 


^(^o)  =  2^  /  fzT.'^C       •     zo  e  n     z^^v 


(3.5) 


r^/:./w  ■ 


As  in  the  linear-element  CVBEM,  the  contour  integral  can  be  eliminated 
(see  Appendix  B),  and  the  above  equation  reduces  to 


U) 


1      ^' 
^^°^    ^   2^  ^  Cfc  '^*t  +   ^k^-\^k^\   +  Cfc  +  2  <^t  +  ; 
j  =  1 

A:  =  2j  -  1 


W 


here 


(~fc  +  2^fc-n)(~fc  +  2  ~  ~fc  +  l)^fc 


D 


(3.6a) 


(^fc  -h  2  +  ^fc  +  l)(^fc  -t-  2  -  ^fc  -n)[(^fc  +  2  -  ^fc)  +  ^oM 

D 
D 


) 


(3.6b) 
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Cjt  +  i  — 


D 


D 

/    2  —  z^    ) 

i^k  +  2-^k)y      '"'2  +^o[(^t  +  2-^fc)  +  ^0^ 

D 


) 


(3.6c) 


Cjt  +  2     — 


D. 


(-fc  + 1  +  ^kK^k  + 1  -  2fc)[(^fc  +  2  -  2fc)  +  2o^fc] 


(2fc  + 1  -  -fc)( 


D 

(•^fc  +  2  ~  ^k) 


^  ^oli^k  +  2  -  2fc)  +  ^oh]} 


D 


(3.6d) 


D     =     (2t  +  2-2fc)(^/:  +  l-2fc)(^fc  +  2-2fc  +  l) 


(3.6e) 


/il   =   In 


(~k  +  2  ~  ~o) 


(~fc-^o) 


In 


i^k  +  2~  ^0) 


(^it-^o) 


+  z^(A:  +  2,  il-;  0) 


(3.6f) 


See  Fig.  3-1  for  the  definition  of  e{k  +  2,  it;  0). 

Eqn.  (3.6)  is  analogous  to  Eqn.  (2.6)  in  the  Hnear-element  CVBEM,  and, 
Hke  Eqn.  (2.6),  it  embodies  two  equations — one  for  (f>{zQ)  and  one  for  V'(-2^o)- 
Close  inspection  of  Eqn.  (3.6)  reveals  that  it  may  be  used  to  calculate  the  value 
of  Q  at  any  location  Zq  /  Cj.,  provided  the  values  of  u  are  known  at  all  nodes. 
Thus,  one  can  apply  Eqn.  (3.6)  not  only  at  interior  points,  but  also  at  the 
"middle"  node  of  any  boundary  element,  Tj  {zq  =  2jt  +  i). 

The  quadratic-element  analog  to  Eqn.  (2.7)  which  can  be  applied  at  the 
"end"  nodes  (20  =  ^k)  is  derived  in  Appendix  B,  and  is  given  by 
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1 

2X2 


i^fc  =  0—A  Ffc_2Wfc-2  +  Ffc.iWfc.i  +  Fjto^fc 


+      Ffc  +  l    UJfc  +  1      +      Fjt^; 


2  ^k  +  2 


M  \ 

.   .'.=  1      .  ) 


k  =  2j-l 
I  =  2i~l 


(3.7a) 


where 


^-^   "   I   2(.,_i-.,.2)    J 


(3.7b) 


Ffc-i  = 


(2/t-^fc-2)^ 


l2(^fc-2rjt_i)(^fc_i-Zfc_2)l 


(3.7c) 


F.        =      In 


{^k  +  2-^k) 


{^k-2-^k) 


+  ie{k  +  2,  k-2;  k) 


+ 


3zfc-2zfc_t-z;>._2  _  3zfc-20fc  +  i-Zfc^2l 

2(2fc-2fc-l)  2(2fc-2fc  +  i)  J 


(3.7d) 


F|t  +  i  - 


(^jt  ~  -it +  2) 


[2{Zk  —  2k  +  l){^k  +  l~  ^k  +  2)} 


(3.7e) 


Ffc  +  2  - 


—  Zi^  +  2z^^i  —  ^fc4.2| 
2(Zfc  +  i-.~,  +  2)        J 


(3.7f) 


C/ 


(^/  +  2'^/  +  l)('^f  +  2  ~  ~/  +  l)^i 

D 


(2:/  4. 2  +  2/  + 1 ) (^/  +  2  -  ^/  + 1 ) [(2/  +  2  -  2/)  +  ^fc/i/] 


D 

.2         ^2 


(■2         _     2\ 


D 


(3.7g) 
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C/  +  1   — 


D 


{^1  +  2  +  ^/)(^/  +  2  -  ^/)[(^/  +  2  -  ^<)  +  ^fc^/] 

D 
D 


1 


(3.711) 


C/  +  2     — 


(^f  +  l^/)(^/  +  l  ~  ^l)K 


D 


(^t  +  i  +  ^/)(^/  + 1  -  ^/) [(^/  +  2  -  ^/)  +  Zkhi] 


D 

,2  ,2 


D 


(3.7i) 


D     —     {^l  +  2~  ^l){^l-k-\~  ^l){.^l  +  2-  ^l^-l) 


(3.7j) 


/i/  =  In 


(^<  +  2-^fc) 


111 


(2/ +  2  "•^Jt) 


(2/  -  2fc) 


+  z^(/  +  2, /;  A:) 


(3.7k) 


See  Fig.  3-1  for  the  definition  of  e{k  +  2,  it  -  2;  A;)  and  e{l  +  2,  /;  it). 

As  with  Eqn.  (3.6),  this  complicated  expression  can  be  broken  down  into 
two  equations — one  for  <;6j.  and  one  for  0j.-  This  breakdown  must  be  performed  so 
that  the  explicit,  implicit,  and  hybrid  methods  of  solving  for  the  unknown  values 
of  (f>  and  0  at  the  boundary  nodes  may  be  implemented  just  as  in  the  linear- 
element  CVBEM.    The  desired  nodal  equations  as  derived  in  Appendix  B  are 


h  = 


2-K 


+    FLl<?^fc-l     +    Fj^.i   Vfc-1 
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+  H  +  ih  +  i  +  Fk  +  i^k  +  i 

M 

r=V^/  +    Cf+2<A,  +  2    +    C/V2  0,  +  4       (3.8) 

/  =  2.  -  1  J 


and 


+  F«<?i,  -  F{^,   . 

+  F^+i<l>k  +  i  -  Fi  +  xVfc  +  i 

+    Ffc  +  2  ^A:  +  2    -    Fjt  +  2  V'jt  +  2 


M 

'^=V,-f  ■  +  Cr^2^i  +  2  -  C/^2  0/  +  2]^  (3.9) 

/  =  2i  -  1 


+    Cf+2<^/  +  2    -    C/+2  0/  +  2]|    (3. 


Here,  the  superscripts  R  and  I  refer  to  the  real  and  imaginary  parts  of  the 
superscripted  quantities,  respectively.  Complete  expressions  for  these  real  and 
imaginary  parts  are  given  in  Appendix  B. 

3.1  Boundary  Conditions  and  Solution  Methods 
The  boundary  conditions  and  solution  methods  discussed  in  Sections  2.1 
and  2.2  in  conjunction  with  the  linear-element  CVBEM  are  equally  applicable 
when  quadratic  elements  are  used,  with  some  minor  modifications.  The  four 
types  of  boundary  conditions — Dirichlet,  Neumann,  Robin,  and  stream 
function — and  their  associated  equations  (Eqns.  (2.11)  through  (2.14))  remain 
unchanged  except   for  the  replacement   of  nodal  index  j  with  index  k.      The 
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derivation  of  the  explicit,  implicit,  and  hybrid  method  equations  require  only 
that  Eqns.  (3.S)  and  (3.9)  be  used  in  place  of  Eqns.  (2.8)  and  (2.9)  as  detailed 
below. 

Explicit  Method:   For  a  Dirichlet  condition  imposed  at  node  fc,  one  sets 
4>k  =  't^k  —  4>k  ^-nd  V'fc  =  V'jk  ill  Eqn.  (3.8)  to  obtain 


+  H  +  i<l>k  +  i  +  F^+i^k  +  i 

+    H  +  2<l>k  +  2  +    Fj^+2^fc  +  2 
M 

l'=V.-/  +  G\^2  <f>U2  +  C/V2  V,  +  4       (3.10) 

/  =  2.-l 


This  equation  is  referred  to  as  the  explicit  phi  nodaJ  equation.    Eqn.  (3.9) 
is  dropped. 

For    a    stream    function    condition    imposed    at    node    k,    one    sets 
V'jk  =  ^A  =  V^k  s^^  ^k  =  <f>k  ill  Eqn.  (3.9)  to  obtain 

^k    =      -i|  F^-2<f>k-2    -    H.2^k-2 

+  Fi«_i  <?>,_!  -  Fi_iVjt-i 
+  Ff  ^,  -  F(  ^, 
+  F^+i<Pk  +  i  -  Fi  +  i0fc  +  i 
+  F/,  +  2  ^fc  +  2  -  Ffc  ^  2  0fc  +  2 
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M 


+  E  [cf  <i>i  -  cj  ^i  +  c/Vi  <^,+i  -  cf+i  0,+, 


l'=\-l  +    Cf+2<?^/  +  2    -    Cf+2^/  +  2]f    (3-11) 


/  =  2i  -  1 


This  equation  is  referred  to  as  the  expHcit  psi  nodal  equation.    Eqn.  (3.8)  is 
dropped. 

For  the  explicit  method,  solution  is  effected  by  applying  Eqn.  (3.10) 
at  each  Dirichlet  node  and  Eqn.  (3.11)  at  each  stream  function  node  to 
form  a  system  of  simultaneous  equations  which  is  solved  for  the  unknown 
values  of  (^  and  %J).  . 

Implicit  Method:  For  a  Dirichlet  condition  imposed  at  node  k,  one  sets 
(j)if  =  (f>ic  and  0^  =  ■^jt  in  Eqn.  (3.9)  to  yield 


^k  =    -^1         F«_2<^,_2  -  FLzV-jfe-: 


■2 


+  F«  <f>,  -  F(  v-; 


k 


+  ^k  +  i<l>k  +  i  -  Ffc  +  1  0jt  +  i 

+    Fk+2^k+2    -    Fk+2^k+2 


M 

.  .'  =  1  ] 

fc  =  2ji-l  +    Ci^+2<f>l  +  2    -    QV  2  '/'/ +  2]  r    (3.12) 

l  =  2i-l 


This  equation  is  referred  to  as  the  implicit  psi  nodal  equation.    Eqn.  (3.8) 
is  dropped. 

Similarly,  for  a  stream  function  condition  imposed  at  node  k,  one  sets 

tpic  —  x/^f,  and  <pf.  =  (f)i^  in  Eqn.  (3.8)  to  obtain 
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?fe    =    i|  H-2</>k-2    +    ^k-2^k-2 

+  FLi-^fc-i  +  F^-r^k-1 
+  Fi^,  +  F^  ^, 

+  Fi  +  i<^fc  +  i  +  Fj^+i  Vfc  +  i 

+    H  +  2<l>k  +  2    +    Fk  +  2i^k  +  2 
M 

'^  =Vi  ^  i  +   Cf  +  2  <A,  +  2    +   C/V  2  0/  +  2]|      (3.13) 

This  equation  is  referred  to  as  the  implicit  phi  nodal  equation.    Eqn.  (3.9) 
is  dropped. 

In  a  manner  similar  to  the  explicit  method,  solution  by  the  impUcit 
method  is  effected  by  applying  Eqn.  (3.12)  at  each  Dirichlet  node  and  Eqn. 
(3.13)  at  each  stream  function  node  to  form  a  system  of  simultaneous 
equations  which  is  solved  for  the  unknown  values  of  (f>  and  ip. 

Hybrid  Method:  For  a  Dirichlet  condition  imposed  at  node  k,  one  sets 
(f>k  =  (f>k  on  the  left  side  of  Eqn.  (3.8)  and  <^t  =  <f>,^  and  ip^.  =  '^p^.  on  the  right 
sides  of  Eqns.  (3.8)  and  (3.9)  to  obtain 

^*    ^    i{  H.2<l>k-2    +    'Pk-2^k-2 

+    H-lh-1    +    Fk-l^k-l 

+  Fil,  +  Ff  ^, 

+  H+ih+i  +  F^+i^k+i 

+    H  +  2h  +  2    +    F^+2i>k  +  2 
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+   Y,  [Cj  ct>i  +  Cf  V-/    +  Cf+i  <^,  +  i  +  C/Vi  V'z  +  i 

1;  '^  i  +  cf + 2  <?^/ + 2  +  c/V  2  V'/ + 2]|   (3.14) 


27r 


+  ^^-x4>k-i  -  FLiV'fc-1 
+  F^  ?,  -  Fl  fe 

+  F^+i<^fc  +  i  -  Fi  +  i  Vfc  +  i 

+    Ffc  +  2  ^fc  +  2    ~    Ffc  +  2  V'fc  +  2 


i  =  1 


yv7 


fc  =  2i-l  +    Cj^2  <^/  +  2    ~    C/  +  2V'/  +  2]f    (3.15) 

/  =  2i-l  J 


These  equations  are  referred  to  as  the  Dirichlet  hybrid  phi  nodal  and  psi 
nodal  equations,  respectively,  since  Eqn.  (3.14)  is  explicit  while  Eqn.  (3.15) 
is  implicit. 

For  a  stream  function  condition  imposed  at  node  A:,  one  sets  V'fc  =  V'lt 

on  the  left  side  of  Eqn.  (3.9)  and  V**  =  V'fc  and  (f)^  =  (f>k  on  the  right  sides  of 

•  ."      ■•.      I.    -    >•    ;.  ■  - 

Eqns.  (3.8)  and  (3.9)  to  yield 


2;r 


^k    =    i^  Fi_2<^fc_2    +    F^_2V'ik-2 

+  F[_i<^,_i  +  F^_i^,_i 

+  Fi  ^fc  +  F«  i,, 

+  Fii  +  i<^fc  +  i  +  Fj^+iVjt  +  i 

+  Fi  +  2<^fc  +  2    +    Ff+2  0fc  +  2 


JU 


rS^^  +    C[+2<?^/  +  2    +    C/V2  0/  +  2]}       (3.16) 

/  =  2i  -  1  J 


■i-  's 
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+  F^  ^,  -  Fi  tA, 

+  Fk+ih  +  i  -  F(  +  iVfc  +  i 

+    Ffc  +  2  </>*:  +  2    -    Ffc  +  2  V'fc  +  2 

rS'-i'  +  C/\2  <;i,  +  2  -  C^2 1/.,  +  4  (3.17) 

/  =  2i-l  ^  J 

These  equations  are  referred  to  as  the  stream  function  hybrid  phi  nodal 
and  psi  nodal  equations,  respectively,  since  Eqn.  (3.16)  is  implicit  while 
Eqn.  (3.17)  is  explicit.  ■   ;  ,- 

In  the  hybrid  method,  solution  is  effected  by  applying  Eqns.  (3.14) 
and  (3.15)  at  each  Dirichlet  node  and  Eqns.  (3.16)  and  (3.17)  at  each 
stream  function  node  to  form  a  system  of  simultaneous  equations  which  is 
solved  for  the  unknown  values  of  (j)  and  rp.  .      ;  __    .  '    ' 

.      ■        I  .     ".   .       :•        1.-. 

For  nodes  where  Neumann  or  Robin  conditions  are  specified,  Eqn.  (2.12)  or 
Eqn.  (2.13)  (with  j  replaced  by  k)  will  be  used  in  conjunction  with  the  implicit 
phi  nodal  equation,  Eqn.  (3.13),  in  the  solution  of  unknowns. 

By  applying  Eqns.  (2.12)  and  (3.13)  at  all  Neumann  nodes,  Eqns.  (2.13) 
and  (3.13)  at  all  Robin  nodes,  and  either  the  explicit,  implicit,  or  hybrid  method 
equations  at  all  Dirichlet  and  stream  function  nodes,  a  system  of  linear  algebraic 
equations  of  the  same  form  as  that  for  the  linear-element  CVBEM  is  obtained. 
As  mentioned  previously,  a  detailed  description  of  how  the  systena  is  formed  is 
given  in  Appendix  A.    The  system  can  be  solved  by  any  direct  method  such  as 
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Gaussian  elimination.  Once  the  unknown  values  of  (f)  and  x})  are  found,  one  can 
use  these  values,  together  with  the  specified  values  {(j)  and  i/'),  as  the  boundary- 
nodal  values  needed  by  Eqn.  (3.6)  for  calculating  Q  at  any   interior  point. 

3.2  Additional  Comments 

Here,  a  few  remaining  comments  regarding  the  quadratic  element  CVBEM 
are  made.  As  with  the  linear-element  CVBEM,  the  implicit  method  is 
recommended  over  the  explicit  and  hybrid  methods  for  the  same  reasons 
mentioned  in  Section  2.2.  Also,  the  calculation  of  -J-  at  the  boundary  nodes 
proceeds  just  as  described  in  Section  2.3.  Finally,  corner  nodes  are  still  a  special 
case,  as  will  be  explained  in  Chapter  IV. 

The  description  of  the  quadratic  element  CVBEM  is  now  complete.  In  the 
next  chapter,  a  methodology  for  the  proper  treatment  of  boundary  corners  in 
either  the  linear-element  or  quadratic-element  CVBEM  is  presented. 


\     ( 

■   A 

4 


.  ^iv  ^-  r,  r^Tj" 


CHAPTER  IV 
TREATMENT  OF  CORNERS  IN  THE  CVBEM 


It  was  mentioned  in  Chapter  I  that  the  corners  in  a  domain  boundary  can 
cause  problems  for  the  RVBEMs  since  -^  is  double-vaJued  at  a  corner  node. 
The  CVBEM  avoids  this  problem  by  solving  for  xj)  instead  of  -J-.  Since  ^  is  both 
continuous  and  single- valued  at  a  corner  (just  like  (^),  corner  nodes  can  generally 
be  handled  more  easily  thcin  in  the  RVBEMs.  A  different  "corner"  problem 
arises  in  the  CVBEM,  however,  due  to  the  way  that  Neumann  and  Robin 
boundaxy  conditions  are  implemented.  In  fact,  this  problem  is  not  limited  to 
corner  nodes  but  exists  at  nodes  which  lie  at  any  interface  between  two  different 
boundary  conditions  as  well. 

This  chapter  is  concerned  with  the  treatment  of  corner  nodes  in  the 
CVBEM,  but  due  to  the  nature  of  the  method,  it  is  instructive  to  first  consider 
how  nodes  which  lie  at  boundary  condition  interfaces  (BCI  nodes)  should  be 
handled.    The  principles  developed  for  dealing  with  BCI  nodes  axe  then  used  to 

establish  a  methodology  for  the  treatment  of  corner  nodes. 

•  -•■■■*.■.. 

4.1  BCI  Nodes 
Before  considering  how  BCI  nodes  should  be  treated  in  the  CVBEM,  it  is 
necessaxy  to  clearly  define  what  a  BCI  node  is.  A  BCI  node  is  defined  as  the 
first  node  imposed  with  a  different  type  of  boundaxy  condition  than  the  node  or 
stretch  of  nodes  immediately  preceding  it.  Here,  the  direction  of  travel  is 
assumed  to  be  the  same  as  that  of  the  contour  integral  in  Eqn.  (2.2).  Each  time 
the  boundaxy  condition   changes,   a  BCI  node  is  encountered.      Thus,  in  the 
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sample  domain  of  Fig.  4-1,  there  are  four  BCI  nodes  along  the  boundary.  (The 
letters  D,  N,  and  R  are  used  to  represent  Dirichlet,  Neumann,  and  Robin 
conditions,  respectively,  imposed  at  the  nodes.)  In  this  way,  the  "new" 
boundary  condition  is  assumed  to  begin  along  the  boundary  just  after  the  node 
preceding  the  BCI  node.  Thus,  only  one  condition  is  actually  imposed  at  the 
BCI  node.  A  "doubly-specified"  BCI  node  where  two  boundary  conditions  are 
available  for  specification  is  a  different  case  and  will  be  discussed  in  the  next 
section  along  with  corner  nodes. 

Proper  treatment  of  BCI  nodes  in  the  CVBEM  requires  examination  of  the 
equations  for  implementing  the  boundary  conditions.  Dirichlet  BCI  nodes  pose 
no  problem  whatsoever  as  the  potential  is  directly  specified  via  Eqn.  (2.11).  In 
contrast,  the  equations  for  implementing  Neumann  or  Robin  boundary  conditions 
(Eqns.  (2.12)  and  (2.13))  must  be  applied  with  caution  since  they  are  "upstream" 
equations,  meaning  that  they  require  information  from  the  previous  "upstream" 
node.  For  example,  suppose  that  node  m  is  the  BCI  node  which  starts  a 
Neumann  or  Robin  stretch  of  boundary.  To  properly  apply  Eqn.  (2.12)  or  (2.13), 
the  values  of  [-^)  or  [H.{4>^  —  (f))]  must  be  known  at  node  m  —  1.  If  node  m  —  1  is 
a  Neumann  or  Robin  node,  then  these  values  will  be  known,  and  BCI  node  m 
can  be  treated  like  any  other  Neumann  or  Robin  node;  however,  if  node  m  —  1  is 
a  Dirichlet  node,  then  these  values  will  not  be  known,  and  neither  Eqn.  (2.12) 
nor  Eqn.  (2.13)  can  be  applied  directly  at  node  m. 

One  way  to  remedy  this  deficiency  is  to  treat  node  m  as  the  psi  reference 
node  mentioned  in  Section  2.4.  Equation  (2.12)  or  (2.13)  would  then  be  applied 
at  all  of  the  "downstream"  nodes  along  the  remainder  of  the  Neumann  or  Robin 
stretch  of  boundary  (nodes  m  +  1,  m  +  2,  etc.).  Unfortunately,  only  one  psi 
reference  node  is  allowed  along  the  boundary  since  the  constant  of  integration 
can  only  be  specified  once.    Consider  the  case  shown  in  Fig.  4-2.    Node  3  is  the 
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-     »        ^ 


FIGURE  4-1.   Example  showing  the  location  of  the  BCI  nodes. 
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is  not  known  here 


second  Neumann  BCI  node 


,  ^    ^  psi  reference  node 
'^  '      (first  Neumann  BCI  node) 


V  s 


FIGURE  4-2.      Example    showing    a   section   of  boundary    containing    two 
Neumamn  BCI  nodes. 
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BCI  node  which  starts  the  first  Neumann  stretch  of  boundary,  so  it  becomes  the 
psi  reference  node.  Node  8,  which  is  the  BCI  node  starting  the  second  Neumann 
stretch  of  boundaxy,  cannot  then  also  be  made  the  psi  reference  node  because  the 
reference  value  for  0  hcis  already  been  specified  at  node  3.  Neither  can  the 
Neumann  equation  for  psi  (Eqn.  (2.12))  be  applied  at  node  8  since  ^  at  node  7 
is  not  known. 

There  are  two  possible  remedies  for  the  situation:  (1)  both  -^  and  <f)  must 
be  specified  at  node  7  or  (2)  ^  must  be  estimated  at  node  7  using  the  Cauchy- 
Riemann  conditions  and  a  finite-difference  expression,  e.g., 

(d^  _  (d^jh    ^     07-06 
\dnP        ^dsP         I  z,  -  zg  I 

-     )  " 

This  leads  to  the  following  altered  version  of  Eqn.  (2.12)  which  can  be  applied  in 
lieu  of  that  equation  at  node  8:  ■  -  - 

In  a  more  general  form,  for  node  j,  the  above  equation  becomes 

Other  more  accurate  finite-difference  expressions  can  certainly  be  used  to 
approximate  f^ji  - 1  when  possible. 

If  node  8  had  started  a  Robin  stretch  of  boundaxy,  a  similar  situation 
would  occur,  and  the  two  possible  remedies  would  be  that  (1)  both  Yicj)^  and  (f) 
must  be  specified  at  node  7  or  (2)  H((^  —  <j)^  must  be  estimated  at  node  7  using 
the   Robin   equation,   the   Cauchy-Riemann   conditions,   and   a  finite- difference 
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expression,  e.g., 


This  leads  to  the  following  altered  version  of  Eqn.  (2.13)  which  can  be  applied  in 
lieu  of  that  equation  at  node  8: 

V's  =  V'7  +  \{^{<i>-<i>oo)s  +  it^I^^^i  ]\z,-zA 

In  a  more  general  form,  for  node  j,  the  above  equation  becomes 

Once  again,  other  more  accurate  finite-difference  expressions  can  be  used  to 
approximate  H((^  —  (?ioo)j  _  1  where  appropriate. 

There  is  yet  another  reason  why  a  BCI  node  which  begins  a  Neumann  or 
Robin  stretch  of  boundary  should  be  treated  a^s  the  psi  reference  node  when 
possible.  Consider  a  case  where  a  Neumann  or  Robin  condition  exists  along  some 
paxt  of  the  boundary,  say  from  node  j  to  node  j  -\-  A;;  see  Fig.  4-3.  Equations 
(2.12)  and  (2.13)  show  that,  in  either  case,  the  value  of  V'  at  nodes  j  +  1  through 
j  -\-k  depend  intimately  upon  the  value  of  xj^j.  In  fact,  any  error  in  xj^j  will  be 
carried  over  as  errors  in  0^  + 1  through  ij^j  ^  fc.  The  errors  in  0  along  the 
Neumann  or  Robin  stretch  of  boundary  will,  in  turn,  cause  errors  in  the 
estimated  values  of  cj)  and  V'  at  all  of  other  nodes  along  the  boundary  due  to  the 
coupled  nature  of  the  equations  in  matrix  equation  (2.23).  If  the  Neumann  or 
Robin  part  of  the  boundary  comprises  a  significant  portion  of  the  total  boundary 
length,  this  can  result  in  a  large  amount  of  error.    Thus,  it  is  imperative  that 
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Neumann  or  Robin  condition 
imposed  here 


'BCI  node 


»    ii^*    J         ; 


4.  t     .' 


:i 


FIGURE  4-3.    A  Neumann  or  Robin  stretch  of  boundary. 
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the  value  of  ip  at  the  first  node  of  a  Neumann  or  Robin  stretch  of  boundary  be  as 
accurate  as  possible. 

With  this  in  mind,  the  following  guidelines  are  provided  for  the  choice  of 
the  psi  reference  noded  location: 

1)  If  portions  of  the  boundary  (or  the  entire  boundary)  are  imposed  with 
Neumann  or  Robin  conditions,  let  the  psi  reference  node  be  the  BCI  node  of  the 
longest  Neumann  or  Robin  stretch  following  the  direction  of  contour  integration 
of  Eqn.  (2.2). 

2)  If  the  entire  boundary  is  imposed  with  a  Dirichlet  condition,  no  BCI   f 
nodes    exist,    and    the    optimum    location    of   the    psi    reference   node    can   be 
determined  by  trial  and  error.    Experience  has  shown  that  its  placement  in  this 
"all  Dirichlet"   case  makes  little  difference  in  the  final  overall  solution.     The    ■ 
position  is  not  crucial  to  accuracy.  ■  " 

To  summarize,  Dirichlet  BCI  nodes  require  no  special  treatment  and  can 
be  handled  just  like  any  other  Dirichlet  node.  On  the  other  hand,  Neumann 
and  Robin  BCI  nodes  should  be  treated  using  the  remedies  and  guidelines 
detailed  above,  subject  to  the  limitation  that  only  one  psi  reference  node  is    _. 

allowed.  "  •    ■  '  N 

'.  '  .( 

4.2  Corner  Nodes  And  Doubly-Specified  BCI  Nodes 
When  using  the  CVBEM  to  solve  potential  problems  with  corners  in  the 
boundary,  the  same  nine  traditional  boundary-condition  combinations  in  Table 
1-1  are  possible.  As  in  Chapter  I,  the  optimal  case  where  both  the  upstream  and  ' 
downstream  conditions  are  known  at  the  corner  is  considered.  Non-optimaJ-case 
corner  nodes  (i.e.,  corner  nodes  which  are  imposed  with  only  one  boundary 
condition)  can  be  treated  in  the  same  manner  as  the  BCI  nodes  described  in  the 
preceding  section;  however,  corner  nodes  where  two  boundary  conditions  are 
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imposed  simultaneously  can  be  treated  specially  in  the  CVBEM.  For  such 
"doubly-specified"  nodes,  the  unknowns  in  Cases  1  through  9  are  listed  in  Table 
4-1. 

Since  the  CVBEM  "corner"  problem  is  due  to  the  way  that  Neumann  and 
Robin  boundary  conditions  are  implemented  and  is  not  due  to  the  double-valued 
nature  of  the  normal  gradient  of  the  potential  (as  in  the  RVBEM),  "doubly- 
specified"  BCI  nodes  (i.e.,  BCI  nodes  where  two  boundary  conditions  are 
available  for  specification)  may  also  be  treated  as  "corner  nodes".  From  this 
point  on,  no  distinction  will  be  made  between  "doubly-specified"  corner  nodes 
and  "doubly-specified"  BCI  nodes,  and  the  term  corner  node  will  be  understood 
to  encompass  both. 

Earlier  in  Section  2.1,  four  boundary  conditions  were  discussed,  namely, 
Dirichlet,  Neumann,  Robin,  and  streajn  function.  A  nomenclature  is  now 
introduced  whereby  a  paxticular  node  is  identified  by  the  boundary  condition 
imposed  there.  Thus,  four  types  of  nodes  are  identified — Dirichlet  nodes, 
Neumann  nodes,  Robin  nodes,  and  stream  function  nodes.  A  fifth  type  of  node — 
the  completely  specified  node — was  discussed  previously  in  Section  2.3.  Now,  a 
question  arises  as  to  whether  in  each  of  the  nine  different  corner  cases,  the  corner 
node  can  be  treated  as  one  of  the  five  nodal  types  above.  In  order  to  address  this 
question,  the  equations  associated  with  the  different  boundary  conditions  must  ■ 
again  be  examined  as  they  were  for  BCI  nodes.  The  application  of  these 
equations  at  cirbitrary  corner  node  m  will  now  be  investigated. 

As  mentioned  in  the  Section  4.1,  the  Neumann  and  Robin  equations  (Eqns. 
(2.12)  and  (2.13))  are  "upstream"  equations  for  calculating  tp^,  meaning  that 
they  require  information  from  the  previous  upstream  node,  m  — 1.  Due  to  this 
"upstream"  nature,  applying  either  Eqn.  (2.12)  or  (2.13)  at  a  "doubly-specified" 
comer  node  which  is  the  first  node  of  a  Neumann  or  Robin  stretch  of  boundary  is 
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TABLE   4-1.     The    unknowns    and    corresponding  number    of    equations 

available   in    the   implicit    CVBEM  for    the   nine   possible 

"doubly-specified"       corner       node  boundary       condition 
combinations  listed  in  Table  1-1. 


Number  of 
Case  Unknowns  Eqns.  Available 


1 

,  V' 

2 

2 

^ 

2 

3 

(f>  and  V' 

.      3 

4 

-      -^'    ,. 

2 

5 

^ 

2 

6 

(j)  and  V" 

3 

7 

(f)  and  tp 

3 

8 

(f)  and  0 

3 

9 

^ 

1 
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undesirable  since  it  forces  the  downstream  Neumann  or  Robin  condition  back  to 
the  previous  node  m  —  1  where  the  condition  is  not  imposed  (see  Fig.  4-4). 
Notice  that  this  situation  is  different  from  that  described  in  the  previous  section 
for  BCI  nodes.  The  BCI  node  boundaxy  condition  was  actually  assumed  to  begin 
just  after  the  node  immediately  preceding  the  BCI  node.  Here,  for  the  "doubly- 
specified"  case,  the  upstream  and  downstream  boundary  conditions  are  assumed 
to  terminate  and  begin,  respectively,  at  the  corner  node  itself  (see  Fig.  4-5). 
Thus,  in  Cases  1,  3,  4,  6,  7,  and  8,  although  two  conditions  are  available,  the 
corner  node  should  be  treated  using  the  upstreaxa  condition  rather  than  the 
downstream  Neumann  or  Robin  condition. 

In  contrast  to  the  Neumann  and  Robin  equations,  the  equation  for 
imposing  a  Dirichlet  boundary  condition  (Eqn.  (2.11))  involves  direct 
specification  of  <f>^.    Therefore,  Case  9  is  easily  handled  as  a  Dirichlet  node. 

This  leaves  Cases  2  and  5  which  axe  grouped  together  as  an  upstream 
equation  for  V'm  followed  by  a  direct  specification  of  (f>^.  For  these  cases,  (f)^  can 
be  specified  directly  using  Eqn  (2.11),  and  the  upstreaxn  Neumann  or  Robin 
equation  (Eqn  (2.12)  or  (2.13))  can  then  be  applied  at  node  m.  Clearly,  this  is 
not  one  of  the  five  nodal  types  described  above.  This  new  type  of  node  will  be 
referred  to  as  a  Neumann-Robin/Dirichlet  node.  " 

To  summarize,  it  has  been  shown  that,  of  the  nine  corner  cases,  seven  can 
be  handled  using  previously  defined  nodal  types.  Only  two — Cases  2  and  5 — 
require  the  addition  of  a  new  type  of  node. 

4.3  Rules  for  the  Nine  Corner  Cases  . 

A  list  is  now  presented  which  gives  the  proper  treatment  for  each  of  the 
nine  cases  in  the  impHcit  CVBEM  at  arbitrary  corner  node  m.  The  equation 
numbers  that  are  given  in  parentheses  refer  to  the  appropriate  nodal  equations 
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downstream  Neumann  or 
Robin  condition 


The  downstream  Neumann 
or  Robin  condition  will 
be  forced  back  to  this 
node  if  Eqn.(  2.12)  or 
(2.13)  is  applied  at 
node  m. 


upstream  boundary 
^     condition 


FIGURE  4-4. 


Diagram   showing   the  effect   of  applying  the  downstream 
Neumann  or  Robin  condition  at  corner  node  m. 


61 


Neumann  condition  begins 
along  boundary  just  after 
this  node  -«. 


BCI  node 


Dirichlet  condition  terminates  at  this  node. 
Neumann  condition  begins  at  this  node. 


FIGURE  4-5.  Diagrams  showing  the  location  where  a  boundary  condition 
change  is  assumed  to  take  place  for  a  BCI  node  and  a 
"doubly-specified"  corner  node. 
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for  both  the  lineax-  and  the  quadratic-element  formulations  of  the  CVBEM. 

Case  1:  Treat  the  corner  node  as  a  Dirichlet  node;  tp^  is  then  the  only 
unknown.  Apply  the  implicit  nodal  equation  for  ip^  (Eqn.  (2.17) 
or  Eqn.  (3.12)). 

Case  2:  Treat  the  corner  node  as  a  Neumann-Robin/Dirichlet  node;  ip^ 
is  then  the  only  unknown.  Apply  the  Neumann  equation  for  V'm 
(Eqn.  (2.12)). 

Case  3:  Treat  the  corner  node  as  an  upstream  Neumann  node.  Both  cf)^ 
and  V'm  ^^^  unknown.  Apply  the  implicit  nodal  equation  for  (^^ 
(Eqn.  (2.18)  or  Eqn.  (3.13))  and  the  upstream  Neumann 
equation  for  tp^  (Eqn.  (2.12)). 

Case  4:  Treat  the  corner  node  as  a  Dirichlet  node;  V'm  is  then  the  only 
unknown.  Apply  the  implicit  nodal  equation  for  0^  (Eqn.  (2.17) 
or  Eqn.  (3.12)).        ,  ,        ,    .  ' 

Case  5:  Treat  the  corner  node  as  a  Neumann-Robin/Dirichlet  node;  xp^ 
is  then  the  only  unknown.  Apply  the  Robin  equation  for  V'm 
(Eqn.  (2.13)). 

Case  Q:  Treat  the  corner  node  as  an  upstream  Robin  node.  Both  <^^  and 
ipm  are  unknown.  Apply  the  implicit  nodal  equation  for  4>^ 
(Eqn.  (2.18)  or  Eqn.  (3.13))  and  the  upstream  Robin  equation 
for  V'm  (Eqn.  (2.13)).    ■ 
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Case  7:  Treat  the  corner  node  as  an  upstream  Neumann  node.  Both  <f>^ 
and  V'm  ^re  unknown.  Apply  the  implicit  nodal  equation  for  <f>^ 
(Eqn.  (2.18)  or  Eqn.  (3.13))  and  the  upstream  Neumann 
equation  for  j/'^  (Eqn.  (2.12)). 

Cctse  £:  Treat  the  corner  node  cis  an  upstream  Robin  node.  Both  (j)^  and 
t/}^  are  unknown.  Apply  the  impUcit  nodal  equation  for  (^^ 
(Eqn.  (2.18)  or  Eqn.  (3.13))  and  the  upstream  Robin  equation 
for  V^  (Eqn.  (2.13)). 

Case  9:  Treat  the  corner  node  as  a  Dirichlet  node;  ^^  is  the  only 
unknown.  Apply  the  implicit  nodal  equation  for  xp^  (Eqn.  (2.17) 
or  Eqn.  (3.12)). 

An  exception  to  the  above  nine  rules  involves  the  psi  reference  node 
described  in  Section  2.3.  This  exception  ties  in  with  the  comments  made  in 
Section  4.1  for  BCI  nodes  and  is  listed  below. 

Exception:  The  first  nodal  point  of  the  longest  stretch  of  boundary 
imposed  with  a  Neumann  or  Robin  boundary  condition 
(traveling  from  the  initial  point  along  the  boundary,  keeping 
the  interior  of  the  domain  on  the  left)  should  be  treated  as  the 
psi  reference  node  described  in  Chapter  II. 

This  important  exception  can  supersede  the  rules  listed  for  Cases  1,  3,  4,  6,  7, 
and  8. 

This  completes  the  discussion  of  the  methodology  for  handling  corners  and 
boundary  condition  interfaces  in  the  CVBEM.  Examples  of  solutions  obtained 
using  the  corner  methodology  will  be  given  in  Section  6.1  of  Chapter  VI.    In  the 
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chapter  that  follows,  a  novel  application  of  the  CVBEM  in  the  axea  of  numerical 
grid  generation  is  presented. 


CHAPTER  V 
NUMERICAL  GRID  GENERATION  USING  THE  CVBEM 


As  mentioned  in  Chapter  I,  grid  generation  methods  can  be  divided  into 
two  main  classes — algebraic  methods  and  differential  equation  methods.  The 
algebraic  methods,  being  much  faster,  tend  to  be  the  methods  of  choice  for 
unsteady,  deforming-boundary  problems  where  a  new  grid  is  required  at  each 
time  step.  Unfortunately,  unlike  the  differential  equation  methods,  the  algebraic 
methods  often  do  not  produce  the  smooth,  non-overlapping  grid  lines  essential  to 
obtaining  accurate  FDM  solutions.  ,'  .   .  .--:,■ 

In  this  chapter,  a  variation  of  the  linear-element  CVBEM  which  can  be 
used  to  generate  grids  for  2-d,  simply-connected  spatial  domains  is  presented. 
Since  the  solution  of  Laplace's  equation  is  involved,  this  method  can  be  classified 
as  a  differential  equation  method;  however,  it  is  anticipated  that  the  method 
could  prove  to  be  computationally  efficient  enough  to  bridge  the  gap  between  the 
existing  algebraic  and  differential  equation  methods. 

5.1  Derivation  of  the  Inverted  CVBEM  Equations 
The  linear  element  CVBEM  interior  point  equations  (represented 
collectively  by  Eqn.  (2.6))  will  first  be  rearranged  so  that  if  one  knows  the  values 
of  (f>  and  0  at  an  interior  point,  one  can  solve  for  the  point's  location, 
2o  =  ^0  +  ^Voi  provided  that  the  values  of  4>  and  V"  at  the  boundary  nodes  are 
known.  These  "inverted"  equations  will  form  an  integral  part  of  the  CVBEM 
grid  generation  method.  To  begin,  Eqn.  (2.6)  is  split  into  its  real  and  imaginary 
parts  as 
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N 


^(^o)  =  i  £  {  <^j  +  i[(^o  -  x,)D  +  (yo  -  y,)C] 
+  ^j  +  Aixo-Xj)C  -  {yQ-yj)D] 
-  4>j[(xo-x^^i)B  +  ivo- yj  +  M 


-  ^j[{^o-^j  +  i)C  -  (yo-yj  +  i)^]} 


and 


N 

J  =  l 


^(^o)  =    -  i  £  {        ^j  +  i[(^o  -  ^,)C  -  (t/o  -  2/,)D] 

-  V'i  +  i[(a;o-a^j)D  +  {yo-yj)C] 

-  <l>j[{xo-^j  +  i)C  -  (yo-yj  +  i)D] 
+  ikj[ixo-xj^,)D  +  {yo-yj  +  i)C]} 


where 


C 
D 
F 
A 
B 


=  [M^,  +  i-x,)  +  B{y^^,-y^)]/F 
=  [B(^.  +  i-^i)  -  A(y,  + 1  -  y,)]  /  F 
=  (xj  +  i-x^)^  +  ivj  +  i-yj)^ 


=  In 


i^j  +  i-^o) 


(5.1) 


(5.2) 

(5.3a) 
(5.3b) 
(5.3c) 
(5.3d) 
(5.3e) 


Eqns.  (5.1)  and  (5.2)  can  be  rewritten  as 


N 

E 

i=i 


27r^(2o)=  £{    -^,kC  -  x,  +  iC  -  yoD  +  yj  +  ,D] 

-(l)j[xoD  -  Xj^iD  +  yoC  -  yj +  iC] 

+  0i  +  iKC  -  xf  -  yoD  +  y,D] 

+  <f>j  +  AxoD  -  XjD  +  yoC  -  yyC]} 


and 
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N 

E 

3  =  1 


j  =  1 

-(/>j[xoC  -  :c_,  +  iC  -  ?/oD  +  Vj  +  M 

-^j  +  i[^oD  -  XjD  +  yoC  -  VjC] 

+  <^i  +  ikC  -  xf  -  j/oD  +  yp]  } 

which  can  be  further  rearranged  as 


+  t/o[D(^,-V,  +  i)  +  C(<^,  +  i-<^,)] 

-<^j(-.r^^iD  -  y^  +  iC) 
•   .      +^,  +  i(-x,C  +  yp) 

'■        +4>i^,{-xp  -  y^C)]}  (5.4) 

and 

-2xV^(zo)  =  E{      ^o[D(V',-V',  +  i)  +  C{<f>j^,-<f>j)] 

-^,  +  i(-xp  -  j/^.C) 

+  <^,  +  i(-a:,C  +  yp)]}  -        (5.5) 

Terms  in  the  braces  can  be  expressed  compactly  as  ■ 

Ci,,  =  C(0,  +  i-0,)  +  T){(j)^^,-<t>-)  (5.6a) 

C2,.  =  D(V',-^,  +  i)  +  C(<^,  +  i-<A,-)  (5.6b) 
Ca.i   =    -  0,(-a:,  +  iC  +  y.  +  iD)   -  ,^^.(-x,.^iD  -  y^^.C) 

+  '/',  +  i(-x,C  +  yp)   +  4>,^,{-xp  -  yp)  (5.6c) 
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(5.6d) 


Substituting  Eqns.  (5.6a)  through  (5.6d)  into  Eqns.  (5.4)  and  (5.5)  yields 


N  N 

3  =  1  J  =  1 


N 


i  =  i 


(5.7) 


and 


N  N 

i  =  1  J  =  1 


N 


(5.8) 


With  the  domain  geometry  specified  and  the  values  of  </>  and  rp  known  at 

interior  point  Zq  (the  hats  can  then  be  dropped)  and  at  the  boundary  nodes, 

Eqns.  (5.7)  and  (5.8)  become  two  simultaneous  equations  in  two  unknowns — Xq 

and  j/q.     Unfortunately,  these  equations  are  nonlinear  due  to  the  presence  of 

terms  A  and  B  in  the  summations.     Recall  that  A  and  B  contain  Zq  =  Xq  +  ij/g 

(see  Eqns.  (5.3d)  and  (5.3e)).      ,'•  .  .     •.    ,      ,      i      ..-'   ^   i. 

Since  Eqns.  (5.7)  and  (5.8)  axe  nonlinear,  a  single-point  iterative  solution 

\     ■■'  ■■..'■'■■' 

process  will  be  used  to  solve  the  equations  simultaneously  for  Xq  and  t/q.     To 

facilitate  the  implementation  of  this  process,  Eqns  (5.7)  and  (5.8)  are  combined 

to  yield 


Vo  = 


2n<f>{z,)  -   f;C3,,        27rrP{zo)  -   f;C4,,- 


j  =  i 


N 


+ 


3  =  1 


N 

EC2„ 

3  =  1 


(     N  N 

EC.,.    EC,. 

3=i  ,3=1 

N  '^     N 

ECm    EC2.. 

U  =  1  J  =  1 


(5.9) 
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and 

^0   =   N      '"' ^ Vo  (5.10) 

The  process  of  solving  for  Xq  and  j/q  proceeds  as  follows: 

1)  Guess  initial  values  for  Xq  and  yo- 

2)  Calculate  Cj  j,  Cj,  j,  C3  j,  and  C4  y  using  the  values  of  Xq  and  j/q- 

3)  Calculate  a  new  value  of  tjq  using  Eqn.  (5.9). 

4)  Calculate  a  new  value  of  Xq  using  Eqn.  (5.10);  the  values  of  Cj  _,, 
C2,j,  C3  J,  and  C4  J  are  taken  from  Step  2,  while  the  new  value  of  j/q 
from  Step  3  is  used. 

5)  Repeat  Steps  2  through  5  until  the  percent  difference  between  the  new 
and  previous  values  of  Xq  and  j/g  is  less  than  a  predetermined  constant. 

These  five  steps  comprise  the  desired  iterative  process  for  calculating  the 

location  of  an  interior  point  when  the  values  of  <f>  and  V'  are  known  at  the 

■     ■      f 

boundciry  nodes  and  at  the  point  itself. 

5.2  The  CVBEM  Numerical  Grid  Generation  Algorithm 
The  iterative  process  described  above  is  now  combined  with  the  linear 
element  CVBEM  to  create  a  method  for  generating  grid  systems  in  2-d,  simply- 
connected  spatial  domains.  The  method  will  henceforth  be  referred  to  as  the 
complex  variable  boundary  element  grid  generation  method  or  CVBEGGM  for 
short.  The  fundamental  concept  behind  the  CVBEGGM  is  the  use  of  potential 
lines    and    streamlines    as    grid    lines,    and    it    is    intended    for    domains    whose 


>^•:•i>^^i '  ::\'^. 
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boundaries  can  be  divided  into  four  separate,  continuous,  smooth  or  non-smooth 
^,  ^  ctirves.    Examples  of  such  domains  are  given  in  Fig.  5-1.  , 

The  CVBEGGM  consists  of  the  following  eight  steps: 

Step  1  -  Specify  the  Domain  Geometry.  Prescribe  the  locations  of  the 
nodal  points  along  the  domain  boundary.  The  number  of  nodal  points  will  have 
a  direct  effect  on  the  accuracy  of  the  solution  (i.e.,  the  quality  of  the  grid). 
Generally,  as  more  nodal  points  are  used,  the  grid  quality  improves,  but  the 
amount  of  time  required  to  generate  the  grid  increases  as  well.  There  is  therefore 
a  compromise  to  be  made,  and  a  rule  of  thumb  is  that  one  should  use  cis  maoiy 
nodal  points  as  possible  for  the  computational  time  constraints  imposed. 
Certainly,  the  number  of  nodal  points  should  be  great  enough  to  adequately 
describe  the  boundary  contour  itself  but  need  not  be  equal  to  the  number  of  grid 
points  ultimately  desired  along  the  boundaries. 

Step  2  -  Define  the  Mapping.  The  CVBEGGM  is  intended  to  map  a  2-d, 
simply-connected,  physical  spatial  domain  onto  a  2-d,  rectangular  transformed 
domain.  The  boundary  of  the  physical  domain  will  be  mapped  one-to-one  to  the 

-    .  four  sides  of  the  transformed  domain,  forming  a  so-called  body-fitted  coordinate 

system  [46].  Since  the  grid  lines  of  opposite  families  connect  opposing  sides  of 
the  rectangular  transformed  domain  (see  Fig.  5- 2b),  the  boundary  of  the  physical 

-  domain  needs  to  be  broken  down  into  four  separate  boundary  curves  (see  Fig.  5- 

2a).  Generally,  there  will  be  a  "natural"  way  to  perform  this  brecikup,  but 
consideration  should  be  given  to  the  type  of  grid  desired.  Reference  46  gives  a 
good  discussion  of  various  breakup  considerations  for  simply-connected  domains. 

Step  3  -  Apply  the  Boundary  Conditions.  The  boundary  conditions  which 
should  be  imposed  to  insure  that  the  grid  lines  intersect  the  boundaries 
orthogonally  are  shown  in  Fig.  5-3.     These  boundary  conditions  force  boundary 
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FIGURE  5-1.      Simply-connected  domains  whose  boundaries  can  be  broken 
up  into  four  separate  boundary  curves. 
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b)  The  transformed  domain. 


FIGURE  5-2.    The  physical  and  transformed  domains. 
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FIGURE  5-3.   Boundary  conditions  for  the  CVBEGGM. 
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curves  1  and  3  to  be  streamlines  and  boundary  curves  2  and  4  to  be  potential 
lines.  Streamlines  in  the  domain  will  thus  connect  and  intersect  orthogonally 
boimdary  curves  2  and  4,  while  potential  lines  in  the  domain  will  connect  and 
intersect  orthogonally  boundary  curves  1  and  3.  The  streamlines  and  potential 
lines  will  become  grid  lines  and  will  be  mapped  to  lines  of  constant  ^  and  rf, 
respectively,  in  the  transformed  domain.  These  grid  lines  will  intersect  both 
each  other  and  the  domain  boundaries  orthogonally — excellent  characteristics  for 
a  grid  system  to  possess  [46].  Corner  nodes  should  be  treated  following  the  rules 
set  forth  in  Chapter  IV,  which  are  given  as  follows: 

a)  the  corner  node  between  boundary  curves  4   and   1   is  treated   as  a 
completely  specified  node  ((f)  =  0,  tp  =  0), 

b)  the   corner  node   between   boundary   curves    1    and   2  is   treated   as   a 
completely  specified  node  [(f>  =  l^  rp  =  0), 

c)  the  corner  node  between  boundary  curves  2  and  3  is   treated  as  a 

Dirichlet  node  {(f>  =  1),  and      '  ,,  ■  *  ^ 

."■-/.■  '1  '•   ■ 

d)  the  corner  node  between  curves  3  and  4  is  treated  as  a  Neumann- 
Robin/Dirichlet  node. 

Step  4  -  Solve  for  the  Unknown  Values  of  (p  and  ip  at  the  Boundary  Nodes. 
This  is  accomplished  using  the  linear  element  CVBEM  as  described  in  Chapter 
II.    The  implicit  formulation  is  recommended. 

Step  5  -  Calculate  the  Boundary  Grid  Point  Locations.  There     are     a 

number  of  ways  of  specifying  the  boundary  grid  point  locations,  but  generally, 
one  first  prescribes  equally  spaced  grid  points  along  boundary  curves  1  and  4. 
The  locations  of  these  grid  points  are  calculated  by  approximating  the  arc  length 
along  the  boundary  as  the  cumulative  sum  of  the  lengths  of  the  linear  segments 
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connecting   the   nodal   points.       Linear   spline   interpolation   is    then   used   to 
determine  the  equally-spaced  positions. 

In  order  to  facilitate  the  presentation  of  the  grid  generation  process,  some 
new  notation  is  in  order;  see  Figs.  5-4  and  5-5.  Previously  in  the  linear-element 
CVBEM,  the  nodal  points  were  numbered  from  1  to  N,  proceeding  around  the 
boundary  of  the  domain  from  an  arbitrary  staxting  position  in  the  direction  of 
the  contour  integral  of  Eqn.  (2.2).  Now,  the  nodal  points  associated  with  each  of 
the  four  boundary  curves  will  be  numbered  from  1  to  N„,  with  N^  representing 
the  number  of  nodal  points  along  boundary  curve  m  (m  being  either  1,  2,  3  or  4). 
The  location  of  nodal  point  k  along  boundary  curve  m  will  be  identified  as 
(X^jt^Ym.fc);  see  the  example  in  Fig.  5-4.  Values  of  the  potential  and  stream 
function  at  this  nodal  point  will  be  designated  as  <f>m,k  ^^'^  '4^m,ki  respectively. 
Approximate  arc  lengths  along  the  boundary  at  nodal  point  k  along  boundary 
curve  m  will  be  represented  by  S^  jt)  ^^^  reference  point  for  zero  length  will  be 
(X^  i,Y„,  i).  Note  that  k  takes  on  values  between  1  and  N„  on  each  boundary 
curve  and  that  each  corner  point  is  associated  with  two  boundary  curves. 

In  addition  to  the  nodal  points,  there  will  be  grid  points  located  along  each 
of  the  boundary  curves.  In  general,  these  grid  points  will  not  coincide  with  the 
nodal  points  except  at  the  corners  and  will  always  be  referenced  separately.  The 
number  of  grid  points  along  streamlines  will  be  designated  as  IL,  while  the 
number  of  grid  points  along  potential  lines  will  be  designated  as  JL.  Grid  point 
locations  will  be  given  by  \x{i,j)^  y{iij))  where  i  runs  from  1  to  IL  and  j  runs 
from  1  to  JL;  see  the  example  in  the  physical  domain  shown  in  Fig.  5-5.  Values 
of  the  potential  and  the  stream  function  at  such  grid  points  will  be  designated  as 
(f>{i,  j)  and  0(i,  j),  respectively.  Finally,  approximate  arc  lengths  along  the 
boundary  at  the  grid  points  will  be  identified  as  s{i,j)  where  i  and  _;'  will 
correspond   to   the   i   and   j   indices   of  the   point's   ordered   pair   designation, 
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FIGURE  5-4.   The  notation  used  for  the  nodal  points  in  the  CVBEGGM. 
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FIGURE  5-5.   The  notation  used  for  the  grid  points  in  the  CVBEGGM. 
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(x{i,j),  y{i,j)\  the  reference  point  for  zero  axe  length  along  each  boundciry  curve 
will  be  given  later. 

Using  the  variables  defined  above,  the  locations  of  the  grid  points  along 
boundary  curve  1    are  given  by 


x(z,  1)   =  ^-^ c^(Xi,/  +  i-Xi  ,)  +  Xi,  (5.11) 

and 

yii,  1)   =  f'^^~c^-'(Yi./-fi-Yi./)  +  ^u  (5.12) 


where  i  =  1,  2,  ...,  IL.    Similarly,  the  locations  of  the  grid  points  along  boundary 
curve  4  are  given  by 


a;(i,  i)  =  c— ^ — rQ^(^4,n  +  i-X4,„)  +  X4,„  (5.13) 

'-'4,n  +  l        '-'4,n 

and 

.  ■  >' 

j/(l,;)   =  iiilil^|^(Y4.„  +  i-Y4.J  +  Y4.„  (5.14) 

'-'4,n  +  l        '-'4,n 


where  i  =  1,  2,  ...,  JL.  Here,  indices  /  and  n  identify  the  nodes  along  boundary 
curves  1  and  4  that  immediately  precede  boundary  grid  points  (x{i,l),  y(«,  1)]  and 
(x{l,j),  y{l,j)\  respectively.  The  values  of  /  and  n  are  found  via  a  simple 
ordered  table  search  of  the  approximate  arc  lengths  along  the  boundary  curves  at 
the  grid  points  and  nodal  points.  The  approximate  arc  lengths  for  the  nodal 
points  are  calculated  by  setting 

S„.,i=  0  (5.15a) 

Then 

S.,.=   i:^(X„,.-X„,._i)2  +  (Y„,.-Y^,,_,)2  (5.15b) 
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where  m  =  1,  2,  3,  4,   and  A;  =  2,  3,  ...,  N^.    Similarly,  by  setting 

s(l,l)=0  (5.16a) 

the  approximate  axe  lengths  for  the  remaining  grid  points  along  boundary  curve 
1  can  be  computed  using 


s(i,l)    =    ^^    {xik,  1)  -  x{k  -  l,l)y  +  {y{k,  l)-y{k-l,l))'  (5.16b) 

k  =  2 

and 


s(i-,JL)  =    ^{        {x{k,JL)-x{k-l,JL)Y 

+  {y{k,JL)-y{k-l,JL)yY-'  (5.16c) 


where  i  =  2,  3,  ...,  IL.     Further,  the  approximate  arc  lengths  at  the  remaining 
grid  points  along  boundary  curve  4  can  be   computed  using 


s(l,j)   =    £;>J    {x{l,k)-xil,k-l)y  +  {y{l,k)-y{l,k-l))'  (5.17b) 

k  =  2 

and 

3 


s(IL,i)   =    £{        {x{IUk)  -  x{lL,k  -  I))' 

+  (y(IL,fc)-j/(IL,A:-l))2}o-5  (5.17c) 


where  j  =  2,  3,  ...,  JL. 

Once  the  grid  point  locations  along  boundary  curves  1  and  4  are 
determined,  the  values  of  </>  and  ^  at  the  grid  points  on  these  boundary  curves 
can  be  calculated.  It  was  found  that  linear  interpolation  between  the  values  of  cj) 
and  ^  at  the  boundary  nodes  produced  satisfactory  results  for  (f>  and  ip,  provided 
that  Eqns.  (2.8)  and  (2.9)  were  used  to  calculate  the  boundary  node  values.  The 
equations  used  for  interpolation  axe  given  below. 


so 


■^i.^  +  i  ~^i,i 


(5.18a) 


and 


rkii,!)  =  il!liL_?M(^         _^    )  +  ^  (5.18b) 


where  i  =   1,  2,  ...,  IL  along  boundary  curve  1,  and 

<f>{l,j)     =     f'^^~_l''"ikn^X-kn)     +    Kn  (5.19a) 

'-'4,n  +  l        ^^4,11 

and 


V'(l,  i)     =    q— ^ C^(V'4,„  +  l-V'4,n)    +    V'4,n 

■^4,0  +  1        'J4,n 


(5.19b) 


where  j  =  1,  2,  ...,  JL  along  boundary  curve  4.  The  values  of  /  and  n  have  been 
defined  eajlier,  and  the  approximate  arc  lengths  are  the  same  as  those  calculated 
previously. 

Finally,  the  values  of  <f>  and  xj)  at  the  grid  points  along  boundary  curves  1 
and  4  axe  used  to  calculate  the  grid  point  locations  along  boundary  curves  2  and 
3.  Because  of  the  boundary  conditions  imposed,  each  grid  point  on  boundary 
curve  1  has  a  corresponding  grid  point  on  boundary  curve  3  which  has  the  same 
value  of  <t>.  Likewise,  each  grid  point  on  boundary  curve  4  has  a  corresponding 
grid  point  on  boundary  curve  2  which  has  the  same  value  of  xl).  For  each  grid 
point  on  boundary  curves  1  and  4,  the  corresponding  grid  point  on  boundary 
curves  2  and  3  is  located  by  using  linear  interpolation  between  the  values  at  the 
nodal  points.    The  equations  used  to  perform  this  operation  are  given  below. 


x(i,JL)   =  t^^^^^-^(X3,,^i-X3,,)  +  X3,,  (5.20a) 

and 

y(i,  JL)   =  t^''^^"^^  (Y3,,  +  i  -  Y3,;)  +  Y3.,  (5.20b) 

'P3,/  +  l-<P3,/ 
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where  i  =  1,  2,  ...,  IL  along  boundaxy  curve  3,  and 


x(IL,j)   =  y^-^^     t^''"(X..„^i-X,,„)  +  X2,„  (5.21a) 

and 

t/(IL,  ;•)  =  y'-^^"^"  (Y,,„^j  -  Y,,„)  +  Y,,„  (5.21b) 

^2,n+l-V'2,n 


where  j  =  1,  2,  ...,  JL  along  boundary  curve  2.  This  time,  the  values  of  /  and  n 
are  found  via  a  simple  ordered  table  search  of  the  values  of  <f)  and  ^  at  the  grid 
points  and  nodal  points  along  boundary  curves  3  and  4,  respectively.  Once 
again,  Eqns.  (2.8)  and  (2.9)  should  be  used  to  calculate  the  boundary  node  values 
of  <f)  and  rj^. 

Step  6  -  Calculate  the  Initial  Interior  Grid  Point  Locations.  There  are 
several  ways  which  an  initial  distribution  of  interior  grid  points  can  be  specified; 
here  it  is  recommended  that  transfinite  bilinear  interpolation  be  used  [45].  This 
method  is  quite  computationally  efficient  and  has  been  found  to  provide 
acceptable  initial  guesses  for  the  interior  grid  point  locations.  The  equations  for 
calculating  the  initied  interior  grid  point  locations  using  this  method  are 

x{ij)   =    (1  -r/)  x{i,  1)  +  r;  x{i,  JL)  +  (1  -^  a:(l,  ;)  +  (  x(IL,  j) 
-  [    (l-0(l-7)x(l,  1)+  e(l-'?)x(IL,  1) 
(l-0^x(l,JL)  +  e»/x(IL,JL)        ] 

and 

y{ij)   =    (1-7/)  y{i,  1)   +  Tj  y{i,  JL)  +  (1  -^  2/(1,  J)   +  ^  vi^U  j) 

-[  (i-0(i-r/)y(i,  1)  +  ai-'/)y(iL,  1) 

(l-0'/2/(l,JL)        +  ^riy{IL,JL)         ] 
where  ' 

^        IL-l  '^  ~  JL-1  ,  •  . 
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and 

i  =  2,  3,  ...,  IL-1  j  =  2,  3,  ...  JL-1 

Step  7  -  Calculate  the  Final  Interior  Grid  Point  Locations.  This  is 
accomplished  by  using  the  iterative  procedure  described  in  Section  5.1.  Each 
interior  grid  point  location  is  calculated  using  this  procedure. 

Step  8  -  Calculate  the  Grid  Point  Locations  in  the  Transformed  Domain. 
Having  calculated  the  locations  of  the  grid  points  in  the  physical  domain,  it  is 
now  necessary  to  calculate  their  mapped  locations  in  the  transformed  domain. 
This  operation  could  have  been  performed  earlier,  but  it  is  sufficient  to  do  so 
now.  The  locations  of  the  grid  points  in  the  transformed  domain  are  given  by 
the  ordered  pairs  ((^,,7/j)  where  .  ,  ., 

(i  =  (i-l)  Ae  i  =  1,  2,  ...,IL 

Jlj  =  (i-l)  At/  i  =  1,  2,  ...,  JL 


and 


A^  =  TT^        ■    "  A7,=      1 


IL-1  " ''  ~  JL  -  1 


5.3  Additional  Comments  On  The  CVBEGGM 
This  section  is  included  in  order  to  mention  some  of  the  more  subtle  points 
associated  with  the  CVBEGGM. 

It  was  stated  in  Step  1  of  the  method  that  the  number  of  nodal  points  used 
can  affect  the  quality  of  the  resulting  grid.  If  the  boundaries  are  composed 
mostly  of  flat  segments,  then  very  few  nodal  points  may  be  required  to  describe 
the  domain  geometry.  This  does  not  mean  that  few  nodal,  points  are  required  to 
produce  an  acceptable  grid,  however,  since  the  values  of  (f>  and  0  may  deviate 
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substantially  from  linearity  over  a  geometrically  lineax  boundaxy.  Since  the 
CVBEGGM  is  a  numerical  process,  some  irregularities  in  smoothness  and 
orthogonality  axe  to  be  expected  no  matter  how  many  nodal  points  are  used. 
(One  might  view  the  CVBEGGM  as  an  approximate  numerical  conformai 
mapping.)  If  no  errors  were  present  along  the  boundaxy,  then  only  round-off 
error  would  occur,  but  generally  there  will  be  errors  at  the  boundaxy  due  to  the 
lineax  approximation.  Increasing  the  number  of  nodal  points  usually  increases 
the  accuracy  of  the  solution,  thus  improving  smoothness  and  necir-orthogonality. 
One  word  of  caution,  however.  Adding  nodal  points  also  increases  the  number  of 
simultaneous  equations  which  must  be  solved.  The  recommended  direct  method 
of  solution  (Gaussian  elimination  with  scaled  paxtial  pivoting)  will  be  subject  to 
roundoff  errors  when  the  number  of  equations  (nodes)  becomes  "large."  ("Large" 
will  be  machine  and  language  dependent.)  Thus,  some  care  should  be  exercised 
when  increasing  the  number  of  nodal  points  to  insure  that  the  solution  is  actually 
improved.  .     -•• 

Another  point  worthy  of  note  involves  the  distribution  of  grid  points  in  the 
physical  domain.  Traditional  elliptic  PDE  grid  generation  methods  employ 
Poisson's  equation,  utilizing  the  source  term  to  control  the  distribution  of  grid 
lines  and  grid  points  [46].  Since  the  CVBEM  is  capable  of  solving  Poisson's 
equation  in  certain  cases,  the  use  of  this  approach  in  the  CVBEGGM  is  currently 
being  investigated;  however,  some  measure  of  grid  line  control  can  be  achieved 
by  simply  varying  the  boundary  grid  point  distribution.  Instead  of  prescribing 
equally- spaced  grid  points  along  boundaries  1  and  4  (as  described  in  Step  5),  one 
can  vaxy  the  spacing  using  one- dimensional  stretching  functions.  A  list  of  such 
functions  is  given  in  Reference  43.  Unfortunately,  because  of  the  smoothing 
characteristics  of  the  Laplacian,  grid  lines  tend  to  be  equally  spaced  away  from 
the  boundaxies  regardless  of  the  stretching  function  used  [46].    Still,  vaxiation  of 
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the  boundary  point  distribution  does  provide  a  simple  albeit  limited  way  of 
controlling  the  location  of  the  grid  points  in  the  physical  domain. 

This  concludes  the  description  of  the  CVBEGGM.  Examples  of  grids 
generated  using  this  method  will  be  included  in  the  next  chapter,  along  with 
examples  demonstrating  the  treatment  of  corners  and  the  use  of  quadratic 
elements  in  the  CVBEM. 


.   •    .  '.•  i 
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CHAPTER  VI 
EXAMPLES  AND  RESULTS 


In  this  chapter,  the  techniques  and  methods  developed  in  the  preceding 
three  chapters  are  verified  and  tested.  First,  the  quadratic  element  formulation 
of  the  CVBEM  derived  in  Chapter  III  is  used  to  solve  two  example  problems, 
and  the  results  are  compared  with  those  obtained  using  linear  elements.  Then, 
the  corner  handling  methodology  developed  in  Chapter  IV  is  applied  to  several 
problems  with  varied  geometries.  The  results  obtained  with  the  CVBEM  axe 
compared  to  available  analytical  solutions  and  to  those  of  other  investigators 
using  RVBEMs  where  appropriate.  Finally,  the  CVBEGGM  of  Chapter  V  is 
used  to  generate  grids  for  four  simply-connected  domains  with  irregular 
boundaries.  t  -.     > 

6.1  Quadratic  Element  Examples  and  Results 
The    quadratic    element    formulation    of    the    CVBEM    was    tested    by 
application  to  the  circular  domain  shown  in  Figure  6-1.    Two  separate  complex 
potential  distributions  were  assumed,  namely 


a;  =   z 


2 


(6.1a) 


which  results  in 


(l>  =  x2-2/2  (6.1b) 

^  =  2xy  (6.1c) 
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FIGURE  6-1.    The    circulax    domain    used    to    test    the    quadratic-element 
CVBEM. 
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1^  =  2[x  cos(^)  -  y  sin(^)]  (6.1d) 


and 

u;  =  e"  (6.2a) 

which  results  in 

4>  =  exp(x)  cos(j/)  (6.2b) 

ip  =  exp(x)  sin(j/)  (6.2c) 

■^  =  exp(x)[cos(j/)  cos(^)  —  sm{y)  sin(^)]  (6.2d) 

For  both  distributions,  Dirichlet  boundary  conditions  were  specified  at  ail  32 
boundary  nodes  (including  node  1,  which  was  treated  as  a  completely  specified 
node  with  ij)  taken  to  be  zero).  The  stream  function  and  normal  gradient  of  the 
potential  were  then  calculated  at  each  node  using  the  quadratic-element 
CVBEM.  Also,  (f)  and  V  were  calculated  at  five  internal  points  labeled  "a" 
through  "e"  in  Fig.  6-1.  '       : 

The  distribution  given  by  Eqn.  (6.1)  was  chosen  specifically  for  its 
quadratic  (f)  and  ip.  Therefore,  the  solutions  obtained  using  quadratic  elements 
should  be  exact  for  them  (except  for  round  off  errors).  This  was  indeed  the  case, 
with  the  values  of  <f)  and  ip  at  both  the  boundary  nodes  and  interior  points  being 
accurate  to  at  Iccist  fourteen  decimal  places.  The  results  for  -J^  at  the  boundary, 
while  not  exact,  were  also  good,  with  a  maximum  percent  error  of  0.59  at  nodal 
point  17  and  a  mean  error  of  0.18  percent  over  all  of  the  boundary  nodes. 

The  exponential  distribution  given  by  Eqn.  (6.2)  was  chosen  to  provide  for 
a   comparison   between    the    solutions   obtained   using   the   lineax-element    and 


''^. 
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quadratic-element  formulations  of  the  CVBEM.   The  results  for  the  percent  error 

in  ^  at  the  boundary  nodes  are  given  in  Fig.  6-2.     Here,  it  is  noted  that  the 

percent  error  in  -^  at  the  boundary  is  a  good  indication  of  the  overall  error  in 

the  computation  for  two  reasons.    First,  in  the  CVBEM,  -^  must  be  calculated 

on 

from   the  nodaJ   values   of  xj)  using  finite-difference  formulae.   This  makes  ^ 
generally  less  accurate  thaji  either  (f)  or  0.     Second,  quantities  at  the  boundary 
which  are  not  specified  by  the  boundary  conditions  tend  to  be  less  accurate  than 
^■\^  those  calculated  at  interior  points.    In  fact,  in  the  solution  of  potential  problems 

by  the  CVBEM,  the  maximum  error  in  the  estimated  values  of  (j)  and  ij)  occur  on 
the  boundary  [7]. 

As  shown  in  Fig.  6-2,  the  errors  at  each  of  the  nodal  points  were  reduced 
significantly  by  the  use  of  quadratic  elements,  with  the  maximum  percent  error 
in  ^  being  13.74  for  linear  elements  and  3.79  for  quadratic  elements,  an  almost 
threefold  improvement.  The  results  at  the  interior  points  axe  given  in  Table  6-1, 
and  the  same  improvement  in  the  accuracy  due  to  the  quadratic  elements  is 

apparent.  "  *.  ^     ,    -        ' 

s      ■■  •         :  .    ......4  I     " 

The  computational  times  for  the  lineax-  cind  quadratic-element  solutions 

referred  to  in  Fig.  6-2  were  19  and  35  seconds,  respectively,  on  an  80286/80287 

based   personal   computer.      However,   it   was  found   that   a  quadratic-element 

solution  using  only  16  equally-spaced  nodes  still  possessed  a  smaller  maximum 

error  in  -5—  than  the  32-node  linear-element  solution.    The  run-time  for  this  16- 
an 

node  quadratic-element  solution  was  11  seconds. 

Additional  examples  of  solutions  generated  using  the  quadratic-element 
CVBEM  axe  presented  in  the  next  section  for  domains  with  corners  in  their 
boundaries. 
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FIGURE  6-2.  Percent  error  in  the  normal  gradient  of  the  potential  at  the 
boundary  nodes  for  the  CVBEM  and  RVBEM  solutions  on 
the  domain  shown  in  Fig.  6-1  imposed  with  the  complex 
potential  given  by  Eqn.  (6.2). 
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V,  . 


TABLE  6-1.  The  percent  error  in  (f>  and  r/'  at  the  interior  points  for  the 
CVBEM  solution  on  the  domain  of  Fig.  6-1  imposed  with  the 
complex  potential  of  Eqn.  (6.2). 


Interior 
Point 

Linear  Elements 
%  Error  m  (j)       %  Error  in 

V. 

Quadratic  Elements 
%  Error  in  <?!>       %  Errc«-  in  rl> 

a 

0.00 

-0.98 

0.00 

0.142 

b 

0.18 

-1.23 

-0.002 

0.123     - 

c 

-0.18 

-2.11 

0.005 

0.354 

d 

-3.46 

-1.45 

0.054 

0.216 

e 

3.48 

-0.56 

0.044 

0.075 
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6.2  Corner  Treatment  Examples  and  Results 
In  order  to  assess  the  accuracy  of  the  CVBEM  in  solving  problems  with 
corners  in  the  domain  boundary,  several  examples  were  investigated.  The  ones 
shown  in  Figs.  6-3  and  6-4  were  taken  from  the  paper  by  Walker  and  Fenner  [32] 
referred  to  in  Chapter  II.  For  Corner  Example  1  (Fig.  6-3),  the  potential  and 
stream  function  distributions  over  the  domain  were  taJcen  to  be  •>••■■ 

4>  =  sin(x)cosh(y) 
and 

0  =  cos(a;)sinh(j/) 

while    for    Corner    Example    2    (Fig.    6-4),    the    somewhat    more    complicated 
expressions 

^  =  ln{[(x  -  2)2  +  {y-  2)Y°-'}  +  x' -  y\ 

and     ^  ..,,         ■'•-     .  ^.  : 


i, 


= '^~V^)]  ^ ''^y 


were  assumed.  In  both  examples,  Dirichlet  boundaxy  conditions  were  imposed 
along  all  parts  of  the  boundary  (except  at  one  node,  where  both  <^  and  xj)  were 
specified).  Sixteen  nodes  were  used  in  Corner  Example  1,  while  32  nodes  were 
used  in  Corner  Example  2. 

Corner  Examples  3  and  4  involve  square  domains,  and  they  differ  only  in 
the  boundary  conditions  imposed  as  shown  in  Figs.  6-5  and  6-6.  Here,  the 
potential  and  stream  function  distributions  over  the  domains  were  assumed  to  be 

<!>  =  exp(x)  cos(j/) 
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FIGURE  6-3.    The  geometry  and  nodal  discretization  for  Corner  Example  1. 
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FIGURE  6-4.    The  geometry  and  nodal  discretization  for  Corner  Example  2. 
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FIGURE  6-5.    The    geometry,    nodal    point    distribution,    and    boundary 
conditions  for  Corner  Example  3. 
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FIGURE  6-6.    The    geometry,    nodal    point    distribution,    and    boundary 
conditions  for  Corner  Example  4. 
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and 

V*  =  exp(x)  sin(y), 

and  they  were  used  to  generate  the  conditions  imposed  on  the  boundaries.  These 
examples  were  chosen  to  test  the  CVBEM  corner  methodologies  developed  for 
the  treatment  of  boundary-condition-combination  Cases  1  through  8  in  Table  1- 
1.  (Case  9  has  been  covered  in  Corner  Examples  1  and  2.)  The  four  comers  in 
Corner  Exaimple  3  represent  Cases  1,  3,  5,  and  7,  while  those  in  Corner  Example 
4  cover  Cases  2,  4,  6,  and  8.  Thirty-two  boundary  nodes  were  used  for  both 
examples. 

Finally,  a  2-d  potential  flow  over  a  circular  cylinder  positioned  between 
two  parallel  plates  is  chosen  for  Corner  Example  5  (Fig.  6-7).  This  is  a  more 
prax;tical  problem  which  has  been  solved  using  a  FEM  by  Martin  [55].  The  same 
problem  has  also  been  analyzed  by  Brebbia  and  Dominguez  [28]  and  Alarcon  et 
al.  [30]  using  linear-element  RVBEMs  and  32  nodal  points. 

The  results  for  Corner  Example  1  are  given  in  Fig  6-8.  This  figure  allows  a 
comparison  between  the  errors  in  -ft-  at  the  boundary  nodes  as  calculated  by  the 
linear-element  CVBEM  and  a  quadratic-element  RVBEM.  It  is  noted  that  for 
Corner  Examples  1  and  2,  the  RVBEM  solutions  were  taken  from  the  paper  by 
Walker  and  Fenner  [32],  who  successfully  used  an  "extra  equation  approach"  to 
handle  the  corners  (see  Section  1.1).  As  for  the  present  work,  the  CVBEM 
gradient  values  were  calculated  from  the  stream  function  nodal  values  using  3- 
point  finite-difference  formulae.  At  nodes  2  through  4,  -^  is  zero  so  that  no 
percent  errors  eire  shown. 

As  plotted  in  Fig.  6-8,  the  CVBEM  errors  are  less  than  the  RVBEM  errors 
at  11  of  the  13  boundary  nodes  shown,  and  the  maximum  error  along  the 
boimdary  is  67  percent  lower  in  the  CVBEM  than  in  the  RVBEM.     This  is 
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FIGURE  6-7.    The   problem   description   and  nodal   point   distribution  for 
Corner  Example  5. 
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FIGURE  6-8.  Percent  error  in  the  normal  gradient  of  the  potential  at  the 
boundary  nodes  for  the  CVBEM  and  RVBEM  solutions  of 
Corner  Example  1. 
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encouraging  since  linear  elements  were  used  in  the  CVBEM,  while  quadratic 
elements  were  used  in  the  RVBEM. 

Figure  6-9  presents  the  results  for  Corner  Example  2.  Once  again,  the 
RVBEM  results  are  due  to  Walker  and  Fenner  [32],  but  this  time,  the  CVBEM 
results  were  generated  using  quadratic  elements.  Also,  the  CVBEM  gradient 
values  were  calculated  using  five-point  rather  than  three-point  finite-difference 
formulae.  One  can  see  that  the  percent  errors  in  the  CVBEM  solution  are  less 
than  the  percent  errors  in  the  RVBEM  solution  at  31  of  the  32  nodes,  with  the 
maximum  error  being  0.83  percent  in  the  RVBEM  and  0.32  percent  in  the 
CVBEM.  The  excellent  results  of  both  methods  can  be  attributed  to  the  near- 
quadratic  nature  of  both  the  problem  geometry  and  the  potential  distribution, 
though  the  CVBEM  accuracy  is  clearly  better. 

The  results  for  Corner  Examples  3  and  4  are  presented  in  Figs.  6-10  and  6- 
11.  Along  each  segment  of  the  boundary,  the  percent  error  in  either  (j>  or  -^  is 
shown  graphically.  One  can  see  that  the  results  obtained  using  the  quadratic- 
element  CVBEM  are  good.    In  both  examples,  the  errors  in  the  (j)  were  below  0.2 

O  J.. 

percent,  while  the  errors  in  -^  were  below  4.0  percent. 

In  order  to  illustrate  the  effects  of  changing  the  location  of  the  psi  reference 
node,  Corner  Examples  2  and  3  were  re-worked.  Table  6-2  shows  a  comparison 
of  the  results  obtained  for  Corner  Example  2  with  the  psi-reference-node 
locations  given  in  the  first  column.  It  can  be  seen  that,  in  general,  the  location 
of  the  psi  reference  node  had  little  effect  on  the  results  in  this  all  Dirichlet  case, 
which  is  as  expected.  In  contrast,  Fig.  6-12  shows  the  increase  in  the  error  in  the 
solution  of  Corner  Example  3  when  the  psi  reference  node  was  not  placed  at  the 
beginning  of  the  longest  stretch  of  Neumann  or  Robin  boundary.  The  maximum 
error  in  (f>  increased  from  0.08  percent  to  3.87  percent,  while  the  errors  in  -^ 
increased  significantly  from  2.12  percent  to  159.74  percent. 
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FIGURE  6-9.  Percent  error  in  the  normal  gradient  of  the  potential  at  the 
boundary  nodes  for  the  CVBEM  and  RVBEM  solutions  of 
Corner  Example  2.  >     -■ 
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FIGURE  6-10.  Percent  error  in  the  unspecified  values  of  the  potential  and 
the  normal  gradient  of  the  potential  at  the  boundary  nodes 
for  Corner  Example  3. 
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FIGURE  6-11.  Percent  error  in  the  unspecified  values  of  the  potential  and 
the  normal  gradient  of  the  potential  at  the  boundary  nodes 
for  Corner  Example  4. 


103 


TABLE  6-2.    Results  on   the   boundajy  for   Corner   Example  2  with  four 
different  locations  of  the  psi  reference  point. 


.    d<i> 

Node                                                                                           Comments 
AYig. Max. Avrg. MajL 


Psi  Ref.         %  Error  in  V'  %  Error  in  W^ 


3 

0.00 

0.00 

0.04 

0.32 

arbitrary  node 

13 

0.00 

0.00 

0.04 

0.32 

node  where  -J- ^  0 
on 

17 

0.00 

0.00 

0.04 

0.32 

sharp  corner 

9 

0.00 

0.00 

0.04 

0.31 

non-sharp  corner 

f  .     ;• 
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FIGURE  6-12.    Results  obtained  for  Corner  Example  3  with  two  different 
locations  of  the  psi  reference  node. 
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The  last  results  presented  are  for  Corner  Example  8 — flow  over  a  cylinder 
placed  between  parallel  plates  (see  Fig.  6-13).  No  exact  solution  exists  for  this 
case,  so  the  linear-element  CVBEM  results  are  compared  with  those  of  other 
researchers  using  linear-element  RVBEMs  and  FEMs  in  Table  6-3.  As  one  can 
see,  the  CVBEM  results  compare  well  with  those  of  the  RVBEMs,  with 
differences  less  than  5.4  percent.  Unfortunately,  the  absolute  error  cannot  be 
computed  because  of  the  lack  of  the  exact  solution. 

6.3  CVBEGGM  Examples  and  Results 

The  CVBEGGM  was  tested  by  application  to  the  four  2-d,  simply- 
connected  spatial  domains  shown  in  Fig.  6-14. 

Example  Domain  1  was  chosen  as  a  simple  test  for  the  method  involving  a 
sharp  corner.  The  grid  generated  by  the  CVBEGGM  for  this  domain  is  shown  in 
Fig.  6-15.  The  grid  is  21  x  21,  and  32  nodal  points  were  used  to  define  the 
boundjiry.  One  can  see  that  the  grid  lines  are  smooth  and  that  they  are  nearly- 
orthogonal  to  each  other  and  to  the  boundaries. 

Example  Domain  2  was  chosen  to  show  the  ability  of  the  CVBEGGM  to 
generate  grid  systems  for  domains  with  curved  boundaries.  As  in  Example 
Domain  1,  the  grid  is  21  x  21,  and  32  nodal  points  were  used  on  the  boundary; 
see  Fig.  6-16.  Once  again,  the  grid  lines  are  smooth  and  nearly  orthogonal. 
Notice  that  the  variation  in  spacing  of  the  grid  lines  is  a  natural  consequence  of 
the  method  due  to  the  solution  of  Laplace's  equation  [46].  Grid  lines  are  seen  to 
be  more  closely  spaced  over  convex  boundaries  than  over  concave  boundaries. 

A  more  practical  geometry  is  embodied  by  Example  Domain  3  which 
represents  a  cylinder  in  parallel  flow.  Fifty-five  boundary  nodes  were  used  to 
generate  the  21  x  15  grid  shown  in  Fig.  6-17. 
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FIGURE  6-13.  Results  for  the  the  unspecified  values  of  the  stream  function 
and  the  normal  gradient  of  the  stream  function  at  the 
boundary  nodes  for  Corner  Example  5. 
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TABLE  6-3.    Resiilts    for    ij)    at    nodal    points    12    through    17    of   Corner 
Example  5. 


y-                  tA                     0                 V                     %  Diff- 
Node    coordinate     CVBEM          RVBEM          FEM                (CVBEM  vs. 
(Ref.  3Q) (Ref.  55) RVBEM) 

12  1.00       0.00        0.00      0.00  n/a 

13  1.16     0.4112       0.3898    0.3419        5.34   . 


14 

1.33 

0.7701  ■ 

0.7548 

0.6880 

2.01 

15 

1.54 

1.1751 

1.1667 

1.1072 

0.72 

16 

1.76 

1.5738 

1.5717 

1.5364 

0.14 

17 

2.00 

2.00 

2.00 

2.00 

0.00 
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FIGURE  6-14.   The  four  example  domains  for  the  CVBEGGM. 


109 


t  ^^ 


FIGURE  6-15.   The  CVBEGGM  grid  for  Example  Domain  1. 
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FIGURE  6-16.   The  CVBEGGM  grid  for  Example  Domain  2. 
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FIGURE  6-17.   The  CVBEGGM  grid  for  Example  Domain  3. 
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Another  practical  geometry  is  represented  by  the  final  example  grid— a 
converging- diverging  nozzle  (see  Fig.  6-18).  The  natural  clustering  of  grid  lines 
near  the  nozzle  walls  is  again  a  characteristic  of  the  method.  The  17  x  11  grid 
was  generated  using  52  nodal  points  along  the  boundaxy. 

It  was  mentioned  in  Section  5.3  that  one-dimensional  stretching  functions 
could  be  used  to  control  the  grid-point  distribution  in  the  physical  domain.  An 
example  showing  the  results  of  such  a  procedure  is  shown  in  Fig.  6-19.  The 
geometry  is  that  of  Fig.  6-16,  and  the  grid  lines  have  been  clustered  near  the  left, 
right,  and  upper  boundaries  by  using  stretching  functions  given  in  Reference  43. 

Besides  using  visual  inspection,  the  resulting  grid  systems  were  evaluated 
by  examining  their  metric  distributions.  The  metrics  are  terms  which  axe 
introduced  in  the  transformed  PDE  by  the  mapping  process,  and  their 
smoothness  or  non-smoothness  can  have  a  direct  effect  on  the  quality  of  the 
numerical  solutions  obtained  [43,  44,  46].  Small  slope  discontinuities  can  usually  ' 
be  tolerated,  but  large  ones  can  lead  to  very  poor  solutions.  As  an  example,  the 
metric  distributions  for  Example  Grid  4  are  shown  in  Fig.  6-20.  One  can  see 
that  they  are  nearly-smooth  and  continuous  at  all  grid  points.  These  metrics  . 
were  calculated  nvunerically  from  the  grid  point  locations  using  second-order- 
accurate  finite-difference  formulae. 

As  an  indication  of  the  efficiency  of  the  method,  Example  Grid  3  took  24     " 
seconds  to  generate  on  a  33  MHz  80386/80387  based  persoucd  computer.    Details 
regarding  the  hardware  and  software  utilized  in  this  dissertation  are  given  in  the 
next  section. 

6.4  Hardware  and  Software 
All  of  the  examples  in  this  chapter  were  solved  by  using  either  a  33  MHz 
80386  based  personal  computer  equipped  with  a  33  MHz  80387  math  coprocessor 
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FIGURE  6-18.   The  CVBEGGM  grid  for  Example  Domain  4. 
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FIGURE  6-19.    The     CVBEGGM    grid    for    Example    Domain    2    with 
stretching   functions   used. 
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FIGURE  6-20.    The  metric  distributions  for  Example  Grid  4. 
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or  a  12  MHz  80286  based  personal  computer  fixed  with  am  8  MHz  80287  math 
coprocessor.  The  computer  codes  themselves  were  written  in  Turbo  Pascal 
Version  4.0  in  order  to  effect  an  efficient  implementation  of  the  over  8000  lines  of 
code  on  a  personal  computer.  Three  major  programs  were  written — one  for  the 
linear-element  CVBEM,  one  for  the  quadratic-element  CVBEM,  and  one  for  the 
CVBEGGM — along  with  several  minor  programs  for  data  creation.  Dynamic 
allocation  of  arrays  was  necessary  to  accommodate  the  storage  requirements  of 
the  associated  matrices. 
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CHAPTER  VII 
CONCLUSIONS  AND  RECOMMENDATIONS 


1.  The  CVBEM  for  simply-connected  domains  has  been  extended  to 
include  the  use  of  quadratic  elements  and  interpolating  functions.  This  was 
accomplished  by  introducing  a  second-degree  Lagrange  polynomial  as  the  basis 
fimction  and  then  following  the  derivation  format  used  by  Hromadka  [7]  for 
linear  elements.  The  quadratic-element  CVBEM  nodal-  and  interior-point 
equations  were  found  to  be  quite  similar  in  structure  to  their  linear-element 
counterparts,  albeit  more  complex. 

2.  The  quadratic-element  CVBEM  was  tested  by  application  to  several 
example  problems.  As  expected,  solutions  obtained  using  quadratic  elements 
were  exact  when  a  quadratic  complex  potential  was  imposed  over  the  domain 
(except  for  round-off  error).  Further,  the  use  of  quadratic  elements  in  the 
CVBEM  reduced  the  maximum  error  in  the  solution  of  other  example  problems 
by  at  least  67  percent  over  the  error  associated  with  linear  elements.  This  was 
accomplished  with  an  increased  computational  time  of  approximately  80  percent. 
It  was  also  found  in  one  example  that  the  number  of  nodes  required  for  a 
particular  level  of  accuracy  could  be  reduced  by  50  percent  when  quadratic 
elements  were  used  in  lieu  of  linear  elements.  The  associated  computational 
time  was  11  seconds — a  reduction  of  over  40  percent. 

3.  A  complete  methodology  for  handling  corners  in  the  CVBEM  was 
presented.  It  was  found  that  the  traditional  RVBEM  problem  of  dealing  with 
the  double- valued  normal  gradient  of  the  potential  at  a  corner  does  not  appear  in 
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the  CVBEM  due  to  the  use  of  the  stream  function  tp.  Instead,  problems 
associated  with  the  proper  implementation  of  Neumann  and  Robin  boundary 
conditions  at  a  corner  node  were  encountered  and  resolved  in  a  series  of  nine 
rules  covering  the  nine  possible  combinations  of  Dirichlet,  Neumann,  and  Robin 
conditions  at  a  corner.  One  important  exception  to  the  nine  rules  regarding  the 
placement  of  the  reference  value  of  the  stream  function  was  noted. 

4.  The  methodology  presented  for  handling  corners  using  the  CVBEM  was 
validated  by  application  to  several  example  problems.  It  was  determined  that 
the  CVBEM  was  able  to  generate  solutions  with  greater  accuracy  than  a  well- 
established  RVBEM  [32]  when  boundary  elements  of  the  same  order  were  used. 
This  was  made  possible  by  the  elimination  of  both  numerical  integration  and 
specialized  augmentation  of  the  nodal-equation  matrices  due  to  the  presence  of 
the  corners. 

5.  A  variation  of  the  CVBEM  was  applied  in  the  area  of  numerical  grid 
generation.  This  method,  dubbed  the  complex  variable  boundary  element  grid 
generation  method  (CVBEGGM),  can  be  used  to  generate  grid  systems  for  two- 
dimensional,  simply-connected  spatial  domains.  The  method  consists  of  two 
primary  phases.  First  the  boundary  of  the  domain  is  imposed  with  specific 
boundary  conditions,  and  the  lineax-element  CVBEM  is  used  to  calculate  the 
unknown  values  of  <f>  and  V*  along  the  boundary.  Second,  an  iterative  search  is 
used  to  determine  the  positions  of  the  boundary  and  interior  points  which  are 
located  at  the  intersections  of  designated  streamlines  and  potential  lines.  These 
points  are  then  used  as  grid  points.  The  method  allows  one  to  control  the  grid 
point  distribution  over  two  of  the  four  boundary  curves  and  to  enforce  the  grid 
line  neai-orthogonality  at  both  the  interior  and  the  boundary  of  the  domain. 
The  method  is  similar  to  the  elliptic  PDE  grid  generation  methods  in  that 
Laplace's  equation  is  solved,  but  since  a  boundary  element  method  is  used,  the 
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CVBEGGM  is  more  properly  classified  as  an  integral  equation  grid  generation 
method. 

6.  The  CVBEGGM  was  tested  via  application  to  four  2-d,  four-sided, 
simply-connected  spatial  domains.  It  was  found  that  the  grid  lines  generated 
were  smooth  and  that  they  intersected  both  each  other  and  the  domain 
boundaries  with  near-orthogonaJity.  Grid  spacing  varied  naturally,  according  to 
the  nature  of  the  Laplacian,  with  grid  lines  being  more  closely  spaced  over 
concave  boimdaries  and  more  widely  spaced  over  convex  ones.  The  24  second 
computation  time  for  a  21  x  15  grid  on  a  33  MHz  80386/80387  based  personal 
computer  indicates  that  the  method  is  computationally  efficient. 

Recommendations  for  future  work  include  the  following: 

1.  A  more  detailed  compcirison  between  the  linear-element  and  quadratic- 
element  CVBEM  is  in  order.  Specifically,  one  should  examine  a  wide  variety  of 
problems  to  determine  whether  the  observations  noted  in  this  dissertation  axe 
representative  of  more  general  trends.  It  is  anticipated  that  there  would  be  a 
point  where  the  addition  of  nodal  points  would  begin  to  significantly  increase  the 
computational  time  so  that  the  use  of  quadratic  elements  would  become  more 
attractive.  Further,  as  the  number  of  simultaneous  equations  was  increased,  the 
accuracy  limitations  associated  with  a  direct  solution  method  such  as  Gaussian 
elimination  with  scaled  partial  pivoting  would  be  exposed.  Thus,  there  would  be 
a  practical  limit  to  the  number  of  nodes  allowable  for  acciuracy. 

2.  The  linear-element  CVBEM  has  been  extended  to  the  solution  of 
problems  in  multiply-connected  domains  [21-26].  The  quadratic-element 
CVBEM  should  be  so  extended  as  well. 

4.  The  CVBEGGM  was  derived  using  the  linear-element  CVBEM.  A 
derivation  using  the  quadratic-element  CVBEM  should  be  attempted.  If  such  a 
derivation  proves  successful,  a  study  should  follow  to  compare  the  use  of  linear 
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and  quadratic  elements  in  terms  of  grid  quality  and  computational  efficiency. 

3.  Since  the  CVBEM  has  been  extended  to  the  solution  of  problems  in 
multiply-connected  domains,  it  appears  that  the  CVBEGGM  mjght  be  extended 
to  generate  grid  systems  for  such  domains  also.  Since  this  would  greatly  enhance 
the  applicability  of  the  method,  this  is  a  ta^k  of  great  importance. 

5.  The  CVBEM  is  able  to  solve  the  Poisson  equation  if  a  particular  solution 
is  available.  Since  the  source  term  in  Poisson's  equation  is  often  used  to  provide 
control  over  the  grid  spacing  in  elliptic  PDE  grid  generation  methods,  a  logical 
extension  would  be  to  attempt  this  in  conjunction  with  the  CVBEGGM. 

6.  A  more  detailed  study  of  the  efficiency  of  the  CVBEGGM  should  be 
undertaken.  The  method  should  be  compared  to  existing  algebraic  and  PDE  grid 
generation  methods  to  determine  how  it  stacks  up  both  computationally  and  in 
terms  of  grid  quality.  Computer  programs  used  for  this  purpose  should  be 
written  by  the  same  prograimmer,  coded  in  the  same  language,  compiled  using 
the  same  compiler,  cind  executed  on  the  same  hardware. 

7.  The  grid  systems  generated  by  the  CVBEGGM  in  this  dissertation  were 
evaluated  visually  based  on  grid-line  smoothness  and  near-orthogonality  and  by 
examining  the  smoothness  of  the  metric  distributions.  The  ultimate  test  of  any 
grid  is  whether  accurate  solutions  to  PDEs  can  be  obtained  by  using  it. 
Therefore,  grid  systems  generated  by  the  CVBEGGM  should  be  used  by  finite- 
difference  codes  to  solve  an  array  of  practical  physical  problems.  The  solutions 
obtained  should  be  compared  with  others  which  axe  known  to  be  of  good 
accuracy.  Only  then  can  the  CVBEGGM  be  qualified  as  a  successful  competitor 
in  the  grid  generation  market. 
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APPENDIX  A 
FORMATION  OF  THE  CVBEM  MATRIX  EQUATIONS 

In  this  Appendix,  the  matrix  equations  associated  with  the  expHcit, 
implicit,  and  hybrid  formulations  of  the  CVBEM  are  derived.  The  resulting 
equations  axe  valid  for  both  linear  and  quadratic  elements. 

A.l  Explanation  of  Matrix  Equation  Variables 
The  matrix  equations  derived  in  the  next  three  sections  will  use  a  compact 
notation  for  identifying  the  coefficient  matrices  ajid  the  unknowns.    This  section 
describes  this  compact  notation. 

The  letter  C  will  be  used  to  identify  any  coefficient  matrix.  Two 
subscripts  and  one  superscript  will  distinguish  one  coefficient  matrix  from 
another.  For  instance,  a  particular  coefficient  matrix  could  be  designated  as  C^^. 
The  first  subscript,  j,  is  used  to  identify  the  type  of  equation  that  was  applied  to 
form  the  matrix,  while  the  second  subscript,  k,  identifies  which  nodal  values 
have  the  elements  of  C^^  as  coefficients.  The  superscript,  m,  indicates  the  type 
of  nodal  point  at  which  the  equations  used  to  form  the  matrix  were  applied  and 
is  only  included  when  necessary  for  clarity.  Subscripts  j  and  k  and  superscript  m 
axe  assigned  as  follows: 

j  becomes  e  when  aii  explicit  nodal  equation  is  applied; 
j  becomes  i  when  an  implicit  nodal  equation  is  applied; 
j  becomes  b  when  a  Neumann  or  Robin  equation  is  applied; 
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Thus,  Cj„  represents  the  matrix  of  coefficients  on  the  unknown  nodcd  vcdues  of  4> 
and  V*  which  was  formed  by  applying  an  impUcit  nodal  equation  at  a  series  of 
Neumann  and  Robin  nodes.  The  particular  explicit  nodal  equation  used  and  the 
nodes  at  which  it  was  applied  will  be  clearly  identified  in  the  next  three  sections 
where  necessary.       .       ,  *     ;  .?'*'* 

One  vector  composed  of  constants  which  axe  not  nodal  values  will  appear 
in  the  matrix  equations,  and  it  will  be  designated  as  V.  This  vector  holds  all  of 
the  constant  terms  associated  with  the  Neumann  and  Robin  equations  when  they 
are  applied  at  the  Neumann  and  Robin  nodes. 

Both  the  specified  nodal  quamtities  {(f>  and  ip)  and  the  unknown  nodal 
quantities  {<f>  and  xp)  will  be  represented  as  vectors,  and  the  following  vectors  will 
be  used: 


^  >  ~     specified  values  of  ^  and  tj); 


m 


k  becomes  s  when  the  associated  nodal  value  (<^  or  V*)  is  specified; 

k  becomes  u  when  the  associated  nodal  value  is  unknown  but  its  harmonic 

is  specified;  '"''^, 

k  becomes  n  when  the  associated  nodal  value  is  unknown  and  its  harmonic 

is  also  unknown; 

m  becomes  b  when  an  implicit  nodal  equation  is  applied  at  Neimaann  and 

Robin  nodes; 

m  is  omitted  in  all  other  cases.  (Two  subscripts  eire  sufficient  to  distinguish 

one  matrix  from  another  in  these  cases).  ■ 


a 


{ 
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unknown  values  of  ^  or  V'  at  the  Dirichlet  and  stream  function 
nodes; 

unknown  vaJues  of  (f)  and  ip  at  the  Neumann  and  Robin  nodes; 
all  values  of  <f>„  followed  by  all  values  of  ^„; 

values  oi  (f>  OT  tp  which  have  been  specified  but  which  are  treated 
as  tinknowns  in  the  hybrid  method  formulation. 


The  matrix  formulation  as  presented  is  quite  general  and  assumes  that 
portions  of  the  simply-connected  domain  boundary  have  been  imposed  with 
either  Dirichlet,  Neumann,  Robin,  or  stream  function  boundary  conditions.  In  all 
instances,  it  will  be  assumed  that  there  are  kn  Dirichlet  or  stream  function  nodes 
and  nn  Neumann  or  Robin  nodes  for  a  total  oi  kn  +  nn  =   N  boundary  nodes. 

Equation  numbers  will  be  given  for  the  linear-element  CVBEM  first  and 
then  for  the  quadratic-element  CVBEM.  If  only  one  equation  number  is  given, 
then  the  equation  is  the  same  for  both  linear  and  quadratic  elements. 

A. 2  The  Explicit  Method  Matrix  Equation 
First,  kn  nodal  equations  are  written  using  Eqn.  (2.15)  or  Eqn.  (3.10)  at 
each  Dirichlet  node  and  Eqn.  (2.16)  or  Eqn.  (3.11)  at  each  stream  function  node. 
This  produces  the  following  matrix  equation: 


I  =  ^-Ws  ^  ^' 


0  ^  '"{M 


Then,  nn  nodal  equations  are  written  using  Eqn.  (2.18)  or  Eqn.  (3.13)  at  each 
Neumann  and  Robin  node.    The  resulting  matrix  equation  is 
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{^  =  ci[l]  +  CI 


^.1 


+   C 


la 


(A.2) 


Finally,  nn  Neumann  or  Robin  equations  are  written  using  Eqn.  (2.12)  at  each 
Neumann  node  and  Eqn.  (2.13)  at  each  Robin  node.  This  results  in  the  matrix 
equation 


n 


w  =  c..{|}.c..{g.c..{|j.v 


(A.3) 


Eqns.  (A.2)  cind  (A.3)  can  be  combined  as 
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This  equation  can  be  combined  with  Eqn.  (A.l)  to  yield 
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Equation  (A.4)  represents  the  system  of  equations  which  must  be  solved  for  the 
unknown  values  of  <j)  and  ^  in  the  explicit  formulation  of  the  CVBEM. 
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A. 3  The  Implicit  Method  Matrix  Equation 
First,  kn  nodaJ  equations  axe  written  using  Eqn.  (2.17)  or  Eqn.  (3.12)  at 
each  Dirichlet  node  and  Eqn.  (2.18)  or  Eqn.  (3.13)  at  each  stream  function  node. 
This  produces  the  following  matrix  equation: 


-  'i\ 


{a  -  M^i 


=  c,jry  +  c. 


+  c, 


(S 


(A.5) 


Then,  nn  nodal  equations  axe  written  using  Eqn.  (2.18)  or  Eqn.  (3.13)  at  each 
Neumann  and  Robin  node.   The  resulting  matrix  equation  is 


w  =  ^'-{1}  -  ^-IM  -  '^■"{li 


(A.6) 


Finally,  nn  Neumcinn  or  Robin  equations  are  written  using  Eqn.  (2.12)  at  each 
Neumann  node  and  Eqn.  (2.13)  at  each  Robin  node.  This  results  in  the  matrix 
equation 


{.-„}  =  c..g);c..{|j.c.{|j 


+  V 


Eqns.  (A.6)  and  (A. 7)  can  be  combined  as 
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This  equation  can  be  combined  with  Eqn.  (A.5)  to  yield 
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Equation  (A.8)  represents  the  system  of  equations  which  must  be  solved  for  the 
unknown  values  of  (f>  and  xjji  in  the  implicit  formulation  of  the  CVBEM. 

A.4  The  Hybrid  Method  Matrix  Equation 
First,  kn  explicit  nodal  equations  are  written  using  Eqn.  (2.19)  or  Eqn. 
(3.14)  at  each  Dirichlet  node  and  Eqn.   (2.22)  or  Eqn.   (3.17)  at  each  stream 
fimction  node.   This  produces  the  following  matrix  equation: 


''^ 


=    C, 


:-:}  ^  ^"la 


+  c, 


(A.9) 


Second,  kn  implicit  nodal  equations  are  written  using  Eqn.  (2.20)  or  Eqn.  (3.15) 
at  each  Dirichlet  node  and  Eqn.  (2.21)  or  Eqn.  (3.16)  at  each  stream  function 
node.    This  produces  the  following  matrix  equation: 


S  =  ^"{|j 


^  ^-C-J 


+  c, 


(3 


(A.IO) 


Then,  nn  implicit  nodal  equations  axe  written,  using  Eqn.  (2.18)  or  Eqn.  (3.13) 
at  each  Neumann  and  Robin  node.   The  resulting  matrix  equation  is 


■..It 


w  -  c^lj .  c^g .  0,(1} 


(A.ll) 


127 


Finally,  nn  Neumann  or  Robin  equations  axe  written  using  Eqn.  (2.12)  at  each 
Neumann  node  and  Eqn.  (2.13)  at  each  Robin  node.  This  results  in  the  matrix 
equation 


v^'-^ii 


^^"}  =  ^4l)^^"{lj^^'"{S^^ 


(A.12) 


Eqns.  (A. 9)  and  (A. 10)  can  be  combined  as 
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Eqns.  (A. 11)  and  (A.12)  can  be  combined  as 
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This  equation  can  be  combined  with  Eqn.  (A.13)  to  yield 
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Equation  (A.14)  represents  the  system  of  equations  which  must  be  solved  for  the 
unknown  values  of  4>  ajad  -tp  in  the  hybrid  formulation  of  the  CVBEM. 
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APPENDIX  B 
DERIVATIONS  FOR  THE  QUADRATIC-ELEMENT  CVBEM 

In  this  Appendix,  the  fundamental  equations  presented  in  Chapter  III  for 
the  quadratic-element  formulation  of  the  CVBEM  are  derived.  Section  B.l  deals 
with  the  interior-point  equation,  Eqn.  (3.6),  while  Section  B.2  is  concerned  with 
the  nodai  point  equation,  Eqn.  (3.7),  whose  real  and  imaginary  parts  are  Eqns. 
(3.8)  and  (3.9). 

B.l  Derivation  of  the  Quadratic-Element  CVBEM  Interior- Point  Equation 
In  this  section,  the  steps  required  to  transform  Eqn.  (3.5)  into  Eqn.  (3.6) 
will  be  presented.    By  using  Eqns.  (3.3)  and  (3.4),  the  contour  integral  in  Eqn. 
(3.5)  is  first  rewritten  as  ,  ► 


C-^Q  j  =  l   l(^fc  +  l-2fc)(2fc-2-2fc)  (^fc  +  l-2fc)(^ifc  +  2-2fe  +  l) 

fc  =  2j  -  1 
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(C-^fc)(^fc  +  2-C)  ^  ;■  ,     'Ci 

(^ik  +  1  ~  •^*:)(^fc  +  2  ~-^fc  +  l)  ■  .•: 


By  expanding  out  the  product  terms,  grouping  like  powers  of  (,  and  pulling  all 
constants  out  of  the  resulting  contour  integrals,  Eqn.  (B.l)  becomes 

^2(0    jr    _     ^1  i^k  +  l^k  +  i)  <^k  (2jt-2fe  +  2)  Wfc  +  1 


130 


+ 


i^k^k  +  l)  ^k  +  2 
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The  three  contour  integrals  on  the  right  side  of  Eqn.  (B.2)  can  now  be  evaluated 
analytically  as  follows:  •     .'    "  .    ' 
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By  substituting  Eqns.  (B.3a),  (B.3b),  and  (B.3c)  into  Eqn.  (B.2),  finding  a 
common  denominator,  grouping  the  coefficients  on  Wjt,  ^j^k  +  ii  ^^^  ^k  +  2^  ^^^ 
substituting  the  entire  result  back  into  Eqn.  (3.5),  the  desired  interior-point 
equation,  Eqn.  (3.6),  is  obtained.    It  is  repeated  below  for  completeness. 


U^ 


-M. 
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1       ^ 

1^(20)   =  2^  II  ^fc  ^*  +  ^*  +  i  ^t  +  i  +  Cfc  +  2  '^k  +  : 
i  =  1 

fc  =  2j  -  1 


Ui 


where 


Ct  = 


Cjfe  +  i 


i^k  +  2-^fc  +  l)(^fc  +  2  ~  ■^fc  +  l)^fc 

D 

(gfc  +  2  +  ^fc  +  i)(^fc  +  2  -  ^fc  +  i)[(^fc  +  2  -  ^fc)  +  ^O^fc] 

D 

(^.^2-^.^i){^'''^2~''^  +  "°[(^^  +  ^""^)  +  ^M 

D 


2fc)/»* 


1 


,  (^fc  +  2  +  ^fc)(^fc  +  2  -  ^fc)[(^fc  +  2  -  ^fc)  +  ^o^fc] 


Cit  +  2 


D 


D 


■    -V     .*■ ,      J- 


(3.6a) 


(3.6b) 


(3.6c) 


■•^i 


f^ 


D 


(3.6d) 


D     =    (2fc  +  2-2'ik)(2A:  +  l-2fc)(^fc  +  2-2jt  +  i) 


(3.6e) 


hi,   =  In 


(^Jb  +  2  ~  ^0) 


{^k  -  2o) 


=  In 


(^fc  +  2  "•2^0) 


(^fc  -  ^0) 


+  i9{k  +  2,  k;  0) 


(3.6f) 


-J. 
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B.2  Derivation  of  the  Quadratic-Element  CVBEM  Noded-Point  Equations 
This   section  presents   the   derivation  of  the   quadratic-element   CVBEM 
nodal  equations — Eqns.  (3.7),  (3.8),  and  (3.9).    In  this  effort,  one  takes  the  limit 
of  Eqn.  (3.5)  as  Zq  approaches  Zjt  (see  Fig.  3-1).    To  facilitate  the  limit  process, 
Eqn.  (3.5)  is  first  rewritten  as 


so  that  the  elements  containing  nodal  point  Zj^  (elements  j  —  1  and  j)  may  be 
treated  separately.    Using  Eqns.  (B.2)  and  (B.3),  Eqn.  (B.4)  can  be  rewritten  as 

2mu{zo)   =  T-— w- — r 

\^k       ^k-2)\^k-l        ^k-2) 
{       (2fc-l2fc)M2fc-2o)-ln(^fc-2-^o)] 

-(^*  +  ^ik-l)(^fc-2fc-2) 

-  i^k  +  ^k-  l)2o[ln(2it  -  Zo)  -  M^k  -  2  -  2o)] 

+      ''~2^*'"      +  zoizk  -  Zfc  _  2)  +  zllH^k  -  ^o)  -  H^k  -  2  -  Zo)]} 


{Zk-^k-l)i^k-l-^k-2)  '      ■     '  '      % 

{  -  (^ifc  -  2^fc)H2fe  -  ^o)  -  M^Jt  -  2  - -^o)] 

+  iZk  +  ^k-2){^k-Zk-2)  ■'       '    v;^ 


+  (Zk  +  Zk-  2)Zo[MZk  -  Zq)  -  ln(Zfc  _  2  -  ^o)] 
,2        ,2 


2^        -  zoizk  -Zk-2)-  zlMzk  -  zo)  -  ln(2jt  _  2  -  2o)]} 


■■^y-. 
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-:» 


+ 


(^k 


(^fc-2*:-2)(^fc-2jt-l) 

-{Zk-2  +  ^k-l){^k-^k-2) 

-{Zk-2  +  Zk  -  l)^o[lll(Zfc  -  ^o)  -  H^k  -  2  -  2o)] 

+   ^^'a*"     +  ^o(^fc  -Zk-2)  +  4M^k  -  Zq)  -  ln(2fc - 2  -  ^o)]} 

{       (^fc  +  l^fc  +  2)[ln(2fc  +  2  -  2o)  -  H^k  -  2o)] 

-{zk  +  2  +  Zk  +  i)zo[H^k  +  2  -  ^o)  -  ln(2fc  -  ^o)] 

+  1±±| ^  +  ^„(2:^  +  2  -  2fc)  +  4[ln(2fc  +  2  -  ^o)  -  lll(^ifc  -  ^o)]} 

'     ^    .     ,         ^.  ^fc-n .-,. 

{  -  (^fc^fc  +  2)H2fc  +  2  -  ^o)  -  ln(2fc  -  2o)] 

+  (^fc  +  2  +  ^fc)(^Jt  +  2-2jk) 

+  (2fc  +  2  +  2r*:)2oN2fc  +  2  -  ^o)  -  ln(^fc  -  ^o)] 
/    2         _     2\ 

L!:| L  _  2jj(2^  ^  2  -  2fc)  -  2o[ln(^it  +  2  -  ^o)  -  M^k  -  ^0)]} 

^    ^k  +  2 

{^k  +  2~  ^k){Zk  +  2  ~  •^fc  +  1) 

{       iZkZk  +  x)[Hzk  +  2-^o)-HZk-^o)] 


^■r 


-      m       X 


-(2^fc  +  2fc  +  l)(2fc  +  2-2fc) 
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i^k  +  2ik  +  l)Zo[H2k  +  2  -  2o)  -  H^k  -  ^o)] 

(^fc  +  2-^fc)  ^  ^^(^^  +  2  -  2fc)  +  2o[ln(2fc  +  2  -  2o)  -  ll(^Jk  -  ^o)]} 


•  /i  r 
i#i-i 


(B.5) 


-ft 


Now,  taking  the  limit  as  Zq  approaches  Z/^  in  Eqn.  (B.5)  yields 


t*  .-.-' 

t--^' 


27rzaJfc   =         ^fc-2 


2(2jk-l  —  2jt_2) 


+  <^k-^ 


-{zk-Zk-2f 


2(2fe-2fc-l)(^ik-l-2fc_2) 


+  t^fc 


(,. 


(^fc-2-^fc) 


+ 


3^fc— 2gfc_i — ^fc-; 

2(zfc-2^fc-i) 


3zfc  — 2zfc^.i  —  ■gfc.).2 
2(2jk-2fc  +  i) 


".  *     ^ 


+    Wfc 


+  1 


{Zk-ZkAr2f 


2(2fc-2jt  +  l)(^fc  +  l-'2fc  +  2) 


..^; 


+  '^k 


+  2 


2(2ik  +  l  —  ■2^fc  +  2) 


-tm^ 


— 1 1-> 


(B.6) 


By  utilizing  Eqn.   (3.6)  and  introducing  a  more  compact  notation,  Eqn.  (B.6) 
becomes  the  desired  nodal  equation,  Eqn.  (3.7),  which  is  repeated  here  as 


■■i>« 
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^(^-t)     =    Wi{  F't-2'^fc-2    +    ^k-l'^k-l    +    F 


k^k 


+    Ffc  +  l<^ik  +  l    +    ^k-\-2^k  +  2 


1        ^ 

i  =  1 

i  #  J  -  1 
/  =  2i  -  1 
k  =  2j-l 


(3.7a) 


.  ' 


where 


Fit-2      = 


(3.7b) 


Ffc-i    = 


-(^fc-2fc-2)^ 


2(2fc-2fc-l)(2fc-l-2fe-2) 


(3.7c) 


■f 


=     In 


(^fc  +  2-^fc) 
(^it-2-2fc) 


+ 


3Zfc— 2Zfc_i — Zfc_; 


32fc— 2zfc+i ~ 


2(^fe 


fc  +  l  ~'^fc  +  2  1 
-^it  +  l)  J 


Fjt  +  1    - 


Fjfc  +  2     - 


i^k~  ^k  +  2) 


C, 


-'{ 


_2(2fc-2fc  +  l)(2fc  +  l-^]k  +  2). 

-2fc  +  22fc  +  l-2ik  +  2) 

2(zjt^i -ZH.  +  2) 

(^/  +  2^J  +  l)(^/  +  2-'2^/  +  l)^/ 

(3.7d) 

(3.7e) 

.  (3.7f) 


'^.^■;.-i' 


;^:^y;£, 


■i-i! 


(z/ + 2  +  z\ + 1)(^< + 2  -  ^t + i)[(^i + 2  -  ^t)  +  ZkH 


D 


{^/ 


+  2  -  ^;  +  i){^^^^^V^  +  '^'^f^^'  +  2  -  ^')  +  ^-t/i/]} 


D 


(3.7g) 


C/  +  1  — 
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D 


{Zl  +  2  +  ^l){^l  +  2  -  Zl)[{zi  +  2  -  2/)  +  ^khl] 


(^/  +  2  -  Z,){^^^±^-^  +  Z,[{Z,  ^,-Z,)  +  Z,h,]} 


D 


(3.7h) 


Q+2  - 


D 


D 

f    2        2\ 

(^< + 1  - zi)y' ^\  '''  +  Zkli^i + 2 - ^/)  +  Zkhi]} 

D 


1 


^    —    {Zl  +  2~  Zl){zi  +  1~  ^l)i^l  +  2~  ^t  +  l) 


(3.7i) 


(3.7j) 


vjw? 


/i,  =  In 


{Zl  +  2~  ^k) 
(2/  -  ^k) 


=  In 


(Z,  ^.2  — ^fc) 


(^/  -  ^ik) 


+  idil  +  2,  /;  k) 


(3.7k) 


Equation  (3.7)  must  now  be  split  into  its  real  and  imaginary  parts  to 
produce  Eqns.  (3.8)  and  (3.9).  This  straightforward  albeit  tedious  task  yields  the 
following  results:  ? 


^'^  =  i 


H-2<l>k-2    +    Fi^_2  0fc-2 

+  H-ih-1  +  ^k-ii^k-i 

+  Fi  +  i4>k  +  i  +  F^+i^k  +  1 

+    H^2h  +  2    +    ^k  +  2^k  +  2 


■'/^ 


■% 


-.^^\: 
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.  =  1 

fc  =  2  j  -  1 
l  =  2i-l 


+    C/+2<^/  +  2    +    C/V2^/  +  2]i       (3. 


8) 


and 


^*    =      ~^{  F?-2'^fc-2    -    F(_2Vfc-2 

+  FU,-  Fi  0, 

+  '^k+i<t>k+i  -  H+ii^k+1 

+    Fk  +  2<l>k  +  2    -    Fk^2tpk  +  2 


M 

rS'^'  +  C/V2  <?^/  +  2  -  C^2  0,  +  4  (3.9) 

/  =  2.  - 1  J 


where  '•  >  '      ' 

Ffc^-2  =     (^fc*  (a;fc-i-a^fc-2)-2*  4-1  +  3  *Xfc_i 

*  a^fc-2-4-2  +  (yfc-i-yit-2)  *  (y*:-2  *  yfc-i 

••    .     +  yik-2))  /  (2*(4-i-2*Xfc_i  *Xfc_2    — 

+  4-2  +  (yfc-i-yfc-2)'))  (B.7a) 

H-2  =    -  {xk*iyk-i-yk-2)  +  Xk-i*{yk-2-yk)  +  ^k-2 

*  (j/fc-2/it-i))  /  (2  *  (4-1-2  *Xfc_i  *Xfc_2 

+  4-2  +  (yit-i-yfc-2)'))  (B.7b) 

Ff-i   =     -  (4*(a;fc-i-a;fc-2)  +  4*(-  4-1 

-  Xfc_i  *  a;fc_2  +  2  *  x|_2  +  (yi-i -2/1-2)  *  (j/ib 
+  2/fc-i-2  *yk-2))  +  H  *  (2  *4-i  *  ^k-2 

+  ^k-i  *  iiyk-yk-2)  *  (yfc-4  *  yjfc-i  +  3  *  yk.2) 


-iy 


.  'I 


■M 


m 


-.i?' 
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-  xl_^)-Xk_2*{xl_2  +  yl-2*yk 

*  yfc-i  +  2  *  y^_i-2  *  t/fc-i  *  y*:-2  +  yfc-2)) 
+  4-1  *  {{yk-yk-2? - ^1-2)  +  ^k-i 

*  ^k-2  *  i4-2  +  iyk-2-yk)  *  (3  *  yfc-4  *  j/^.i 
+  yk-2))  +  iyk-yk-i)  *  (4-2  *  (2  *yk-yk-\ 
-  yk-2)  +  iyk-i-yk-2)  *  (j/fc-yfc-2)^))  /  (2 

*{xl-2*  Xk*  Xk_i  +  xl_i+yl-2 
*yik*  2/fc-i  +  J/fc-i)  *  (4-1-2  *  Xfc.i 

*  ^k-2  +  4-2  +  iyk-i-yk-2y)) 


Si 


(B.7c) 


FLi  = 


F?  = 


i4*{yk-\-yk-2)  +  4*ixk-i 

*  (-  yfc-2  *  yfc-i  +  3  *  j/fc_2)  +  a?fc-2  *  (yt-yfc-i)) 

+  ifc*  (2  *xLi  *  {yk-yk-2)-2  ♦Xfc-i  *a:fc_2 

*  (y*-2  *  yi_i+yfc-2)  +  (yfc-2-yjk-i)  *  (4-2 
+  (yfc-2-yfc)  *  (yfc-2  *  yfc_i  +  yfc_2)))-2 

*  xl_i  *  Xfc_2  *  (yfc-yfc-2)  +  ^fc-i  *  (4-2  , 

*  (3  *  yfc-2  *  yk-i-yk-2)-{yk-yk-2?  *  (yk 

-  2  *  j/fc_i  +  j/fc_2))+a:fc_2  *  (yik-i-yfc)  *  (4-2     , 

+  (yfc-2-yfc)  *  (yfc-2*yfc-i  +  yfc_2)))  /  (2 

*{xl-2*Xk*Xk_i  +  xl_i+yl-2 

*yk*yk-i  +  yl-i)  *  {4-1-2  *Xk_i 

*  Xk-2  +  4-2  +  {yk-i-yk-2)^)) 

ln|(2fc  +  2-2lt)  /  (-fc-2-2fc)l 
+    {3*xl  +  Xk*{-    5*Xfc_i-Xfc_2)+2 

*  4-1  +  ^fc-i  *  a;fc-2  +  (yfc-yjt-i)  ♦  (3  *  ^^-2 
*yfc-i-yfc-2))  /  (2*  (4-2*Xfc  *Xfc_i 


(B.7d) 


^ 


..} 
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+  4-i  +  {yk-yk-iY)) 

*  4  +  1  +  ^:^  +  1  *  Xk  +  2  +  iyk-yk  +  i)  *  (3  *  J/fc-2 
*yfc  +  i-yfc  +  2))  /  (2*  (4-2*a;fc*a;fc  +  i 

+  xl+i  +  iyk-yk+iY)) 


(B.7e) 


F{  = 


^(Jk  +  2,  Jt  -  2;  it) 

+  i^k  *  {yk-i-yk-2)  +  Xk-i  *  (yk-2-yk)  +  ^k-2 
*  (yk-yk-i))  /  (2  *  ixl-2  *Xk* Xfc_i 

-  {xk*  {yk+i-yk+2)  +  Xk  +  i  *  {yk+2-yk)  +  ^k+2 
*{yk-yk+i))  /  (2*  ixl-2*Xk*Xk+i 
+  4+i  +  {yk-yk+iY)) 


(B.7f) 


'^k  +  1 


i^k  *  i^k  +  l~  ^k  +  2)  +  ^k  *  i~    ^k  +  1 

-  Xfc^.1  *  Xk+2  +  2  *  xl+2  +  {yk+i-yk+2)  *  (yk 
+  yfc  +  1-2  *  yk+2))  +  Xk  *  (2  *  4+1  *  ^k+2 

+  ^k+i  *  {{yk-yk+2)  *  (yfc-4  *j/fc  +  i  +  3*  j/fc+2) 

-  4+2)-^fc+2  *  (4+2  +  yfc-2*j/fc 

*  yk+i  +  '^  *  yl+i-^  *  Vk  +  i  *  yk+2  +  yl+2)) 
+  4+1  *  {iyk-yk+2y - 4+2)  +  ^k+i 

*  Xk+2  *  (4+2  +  (yjt+2-2/fc)  *  (3  *  yfc-4  *  j/fc+i 
+  yit+2))  +  (2/fc-2/fc+i)  *  (4+2  *  (2  *  Vk-yk+i 

-  yk+2)  +  {yk+i-yk+2)  *  iyk-yk+2)^))  I  (2 

*{xl-2  *Xk*Xk  +  ^+xl^^  +  yl-2 

*  yfc  *  yjfc+i  +  yfc+i)  *  (4+1 -2  *  Xfc  +  i 

*  Xjk  +  2  +  4  +  2  +  (yfe  +  l-J/fc  +  2)^)) 


-^:y-. 


(B.7g) 
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FUx  = 


-  (4*  (yfc+i-yfc+2)  +  4  *  (a:fc+i 

*  {-  yk-^*yk+i  +  ^*yk+2)  +  Xk+2*  {vk-yk+i)) 

4-  Xfc  *  (2  *x|  +  i  *  (yfc-yfc  +  2)-2  *  x^  +  i  *  x^  +  i 

*  (yfc-2  *yfc+i  +  j/jfc+2)  +  (yfc+2-yfc+i)  *  (4+2 
+  (yit+2-j/jt)  *  (yfc-2  *  j/jt+i  +  yit+2)))-2 

*  Xfc  +  1   *  Xfc  +  2  *  (j/fc  ~  2/fc  +  2)  +  ^fc  +  1  *  i^k  +  2 

*  (3  *  yjk-2  *yjt  +  i-yik  +  2)-(yfc-y*:  +  2)^  *  (yt 

-  2  *yfc  +  i  +  yfc  +  2))+Xfc  +  2*  (yfc  +  i-yit)  *  (4  +  2 

+  (yjk+2-yfc)  *  (yfc-2  *  yifc  +  i  +  yfc+2)))  /  (2 

*(xl-2*Xk*  Xk  +  i+xl  +  ^  +  yl-2 

*  yjfe  *  yfc+i  +  yl  +  i)  *  (4+i-2  *  x^+i 
*Xk+2  +  4+2  +  {yk+i-yk+2?)) 


(B.7h) 


^it  +  : 


—  (Xfc  *  (xjt  +  i~^fc  +  2)~2  *  Xjt^.1  +  3  *  Xfc^i 

*  Xfc^.2-x|+2  +  (yfc+i-yfc+2)  *  (yfc-2  *yifc+i 

+    yife  +  2))  /  (2  *  (4  +  1-2  *  Xfc  +  1  *  3:^  +  2  +  4  +  2 

+  (yfc+l-yib+2)^)) 


(B.7i) 


H+2  =   ixk*{yk+i-yk+2)  +  ^k+i*{yk+2-yk)  +  Xk+2 
*  (yk-yk+i))  /  (2  *  (4+1-2  *Xk+i  *  Xk^2 
+  4+2  +  (yfc+i-yifc+2)^)) 


(B.7j) 


Cf  = 


c' 


/   — 


E  = 

.  (NC,)^ 

+  F^ 

.  (NC,)^ 

£2 

+  F' 

E  = 

.  (NC,)^ 

-F* 

(NC,)^ 

E^  +  F^ 


(B.Tk) 


(B.71) 
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'^J  +  i    — 


C/  +  1    — 


pH         - 

^/  +  2    — 


^J  +  2 


E*(NC,  +  i)«  +  F*(NC,  +  i/ 


E^ 

+  F^ 

E^ 

.(NC,  +  i)^ 

-F* 

(NC,  +  i)« 

£2 

+  F2 

E*(NC,^2)'' 

+  F* 

^(NC,  +  2/ 

E^ 

+  F2 

^i 

'(NC,  +  ,)^ 

-F* 

(NC,^^)'' 

2    .    Tr2 


E^  +  F 


(B.7m) 


(B.7n) 


(B.7o) 


(B.Tp) 


(NC,)^  =      -  (2*  A*(x?+,  *(x,  +  2-ijt)  +  X/  +  i  *(-  a:?+2 
+  il  -  2  *  j/,^.1  *  (r/,^2  -  2/fc)  +  y?+2  -  J/fc) 

*y/  +  2-yfc*  (2  *y/+2-yit))  +  a;fc*  (y/  +  i-y/  +  2) 
*  (y/+i  +  y/+2-2  *yfc))-2  *B  *  (x7^.i  *  (y/+2-J/fc) 
+  2  *  x,  +  i  *  (x,  +  2  *  (y/  +  i-y/  +  2)  +  a;fc  *  (y^-yj  +  i)) 

+  a;?+2  *  (yfc-y/  +  i)+2*x,  +  2  *a:fc*  (y,  +  2-J/*) 

+  (y/  +  i-y/+2)  *  (4  +  y/+i  *  (yib-2//+2)  +  y/+2 

*  (x?+i-a:,  +  i  *  3;fc-x?+2  +  x,^.2  *  ifc  + (y/  +  i -y/  +  2) 

*  (y/-J//  +  i-y/  +  2  +  2/fc))-2  *  x,%i  *  a;/  +  2  +  a:;  +  i 

*  (a;?+2  +  2  *  a;/  +  2  *  Xk  +  y]-2  ♦  t/,  *  (2  *  y,  +  i-j/fc) 
+  yi  +  2  *  (4  *  2//  +  i-y/  +  2-2  *  yfc))  +  a;f+2-2  *  a:?+2 

*  Xfc  +  a;,  +  2  *  (-  y?  +  2  *  y/ *  (2  *y,^.2-yfc)  +  2 

*  y?+i-2  *  y/  +  i  *  (y/  +  2  +  yfc)-y/  +  2  *  (3  *  y/  +  2-4 

*  j/jk))  +  2*Xfc*  (y/  +  i-y/  +  2)  *(y/-y/  +  2))  /2  (B.7q) 


.*■ 


(NC,)^  =      -  (2*  A*(x,%i  *(y,  +  2-yfc)+2*x,  +  i  *{xi^, 
*  (y/  +  i-y/  +  2)  +  a;fc*  (yfc-y/  +  i))  +  a:?+2 


>?".* 
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y-M 


*  {yk-yi+i)+'^*  ^1+2*  ^k*  iyi+2-yk)  +  {yi+i 

-  yi+2)  *  {xl  +  vi+i  *  {yk-yi+2)+yk*  iyi+2-yk))) 

+  2  *  B  *  (xf^i  *  ixi  +  2-^k)  +  Xi  +  i  *  (-  xf+2  +  Xk 

-  2  *  y,+i  *  {yi  +  2-yk)  +  yh2-yl)  +  xl+2  *  ^k 
+  ^1+2  *  (-  4-yl+i  +  ^  *yi+i  *yi+2-yk*  (2 

*  y/+2-yfe))  +  ^jk  *  (y/+i-y/+2)  *  (y(+i  +  y/+2-2 

*yi))  +  3;?  *  (y,  +  2-y/  +  i)-2  *  a:,  *  (x,  +  i  *  (i/,-2 

*y/  +  i+yfc)  +  a:/  +  2  *  (- 2//  +  2  *  y/  +  2-yit)  +  a^A: 

*  (y/  +  i-y/  +  2))  +  2  *x?^.i  *  (y,-y,^.2)-2  *!,  +  ! 

*  (x/  +  2*  (2*y/  +  i-y/  +  2-yfc)  +  a:fc*  (y/-y/+2)) 

+  a:?+2  *  (-  2  *  t//  +  y/  +  i  +  3  *  j/;  +  2-2  *  yfc)  +  2 

*x,  +  2  *a:jt  *  (y/  +  y/  +  i-2  *y/+2)  +  (y/  +  i-y/  +  2) 

*  (y? -2  *  y/ *  (y,  +  i  +  y/  +  2-yfc)  +  y/  +  2  *  (2  *y/  +  i    • 

+  y/  +  2-2*yfc)))/2    '    .  ,    •        (B.7r) 


(NC,  +  i)^  =  (2*A*(x?*(x,  +  2-a:fc)  +  a;,  *(-  a;?+2  +  4- 


+  2  +  ^fc~2  ■  /., 

yi*  {yi  +  2~yk)  +  yl+2-yl)  +  xl+2* ^k  +  xi+2 

*  (-  4-y?  +  2  *y/*y,  +  2-yik*  (2  *y/  +  2-yfc)) 

+  Xk*  iyi-yi  +  2)  *  (y/  +  y/  +  2-2  *yjt))-2  *  B  *  (xf    - 

*  (y/+2-yfc)  +  2  *x,  *  (x,+2  *  (y/-y/+2)  +  ^fc*  (y/t  - 

-  y/))  +  a:?+2  *  (yfc-y/)  +  2  *  x,^.2  *  x;t  *  (y/  +  2-yfc)        .  . 
+  (y/-y/+2)  *  (4  +  y/ *  (yik-yj+2)  +  yfc  *  (y/+2 

-  yk)))  +  ^f  +  ^l  *  (-  a;<  +  2-2  *  xjt)  +  x,  *  (-  x^+2         ,•   - 

+  4  *  X,  +  2  *  3;*;  -  3  *  yf  +  2  *  y,  *  (y,  ^.  2  +  2  *  yfc)  .  ; " 

+  y/  +  2  *  (y/  +  2-4  *yfc))  +  a;?+2-2  *  a;?+2  *  a^fc 
+  Xi  +  2  *  (y?  +  2  *y,  *  (y/  +  2-2  *  yk)-yi  +  2  *  (3 

yfc))  +  2*Xfc*(y,-y,  +  2)')/2  (B.7s) 


*  y/+2-4  * 
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(NC,  +  i)^  =   (2*  A  *  (x?  ♦  (y,  +  2-yfc)  +  2  *x,  *  (a;,  +  2  *  (y/-y/  +  2) 

+  ^k*{yk-yi))  +  x'i+2*{yk-yi)  +  '^*^i+2*Xk 

*  (y/+2-yfc)  +  (y/-y;+2)  *  (4  +  y/  *  iyk-yi+2) 
+  yi+2  *  yk-yi))  +  '^  *  B  *  (x]  *  {x,+2-^k)  +  xi 
*{-  a;?+2  +  4-2  *  y/ *  (y/  +  2-yfc)  +  y?+2-yfc)      ^ 
+  xf+2  *  ^k  +  ^i  +  2  *{-  4-2// +  2  *  y,  *  yi  +  2-yk 
*(2*yj+2-yfc))  +  ^fc*  (y/-y/+2)  *(y/  +  y/+2-2 

*yk))  +  xl  *  (3  *y/-y;  +  2-2  *yfc)-2  *x,  *  {xi^2 

*  (yj  +  y/  +  2-2  *  yfc)  +  2  *  x;t  *  (y/-y/  +  2))  +  3;?+2 

*  (-  y/  +  3  *  y;  +  2-2  *  yjt)  +4  *  x,  +  2  *  ^k  *  (yi 

-  yi+2)  +  iyi+2-yi)  *  (y?-2  *y<  *yk-yi+2  *  iyi+2 

-2*y,)))/2'       ."     ■  (B.7t) 

.  .••■'  «■■■ 

(NQ  +  2)^  =    -  (2*  A*(x?*(x,  +  i-Xfc)  +  x,  *(-  x?+i+x^-2         •:. 

•;■     1  *  yi  *  (yj  +  i-yjk)  +  y/  +  i-yD  +  a:?+i  *  a:fc  +  2r,  +  i  .■,, 

"'  *  (-  4-y?  +  2  *  y/ *  y/  +  i-yfc  *  (2  *  y,  +  i-yfc))         ' 

+  arfe*  (y,-y/  +  i)  *(y/  +  y/  +  i-2*yjk))-2*B  *  (x? 

*(y/  +  i-yfc)  +  2  *3;,  *  (x,  +  i  *  (y/-y/  +  i)  +  a;fc*  (y^ 

-  y/))  +  4+i  *  (yfc-y/)  +  2  *x,+i  *Xfc  *  (y,  +  i-yfc) 
+  (y/-y/+i)  *  (4  +  yt  *  (yfe-y/  +  i)  +  y/+i  *yk 

-  yl))  +  ^f  +  ^l  *  (a:/  +  i-2  *  (x,  +  2  +  ^fc))  +  a;, 

*  (  -  2  *  xf  + 1  +  2  *  X,  + 1  *  x;t  +  a:?+  2  +  2  *  X,  ^.  2  *  Xfc 

-  3*yj-2*yi*  (y/  +  i-2  *  y,  +  2-2  *  yfc)4-2 

*  y?+i-2  *  y/  +  i  *  yfc-y/  +  2  *  (y/+2  +  2  *  yfc))  +  2 

*2^/  +  l*^/  +  2  +  ^/  +  l*(—    ^/  +  2  —  2*X,  ^2*^fc  —  y/ 

+  2  *  y,  *  (2  *  y,  +  i-yfc)-j/,^.2  *  (4  *  y,  +  i-y,  +  2 


.^■^: 


^'^v 
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yk)  +  ^k*iyi-yi+2)))  /2  .  (b.7u) 


(NC,  +  2)^   =    -  {2  *  A  *  {xj  *  {y,  +  i-yk)  +  2  *  X,  *  {x,^i  *  {y, 

-  yi+i)  +  Xk  *  {yk-yi))  +  xhi  *  (j/fc-y/)  +  2  ♦  x,+i 
*Xk*  {yi+i-yk)  +  {yi-yi+i)  *  (4  +  j//  *  ivk 

-  J//  +  i)  +  J/fc*  (j//  +  i-J/fc)))  +  2  *  B  *  {xj  *  (x,  +  i 

-  h)  +  xi  *  (-  xhi  +  d-2  *yi*  (yi  +  i-yk) 

+  yl+i-yl)  +  ^l+i  *  Xk  +  ^i  +  i  *  (-  4-y?  +  2 
*yi*  yi  +  i-yk  *  (2  *  yi  +  i-yfc))  +  a:fc  *  {yi-yi  +  i) 

*  (y/  +  y/  +  i-2  *yfc))  +  xf  *  (3  *  y/  +  y,  +  i-2 

*  iyi  +  2  +  yk))  +  2  *x,*  {x,  +  i  *  {y,-2  *  y/  +  i  +  yfc) 
+  a:/  +  2  *  (-  2  *  yi  +  y,  +  2  +  yk)  +  ^k  *  (-  2  *  y, 

+  y/+i  +  y/+2))-2*a;?+i  *(y/-yj+2)  +  2*x,+i 

*  (x,^2  *  (2  *  y/  +  i-y/  +  2-yib)  +  a^fc  *  (y/-y/+2)) 

■-5    +  (2//-y/  +  i)  *  (a:?+2  +  2  *  X/  +  2  *  a;ib-y?-2  *  y, 
i 

•"  *  (y/+i-y;+2-yib)  +  y/+2  *  (2  *y/+i-2//+2-2 

*  J/fc)))  /  2  _       .    , 


!' 


(B.7v) 


E  = 


^/    *  (^/  +  2  ~  ^/  +  l)  +  ^/  *  (^/  +  1  ~  ^/  +  2 

+  (y/  +  i-yj  +  2)  *  (2  *  y/-2//  +  i-2//+2))-a;?+i 

*  a:/  +  2  +  a;,  +  i  *  (x/+2  +  (j// -  ^,  +  2)  *  (^z  -  2 

*y,+i  +  y/+2))  +  a^/+2  *  (y/+i-y/)  *  {yi  +  yi+i 

-    2  *»//  +  2) 


(B.7vv) 


'■4 


F  =  ^?  *  (y/  +  2-y/  +  i)  +  a;/ *  (a^z  +  i  *  (2  *  J//  +  1 

-  2*yi)+  x,  +  2  *  (2*y,-2  *  ^,  +  2)) +  a:?+i 
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*  (y/-j//+2)-2  *  x,  +  i  *  x,+2  *  (j//  +  i-y/+2) 
+  xl+2*iyi+i-yi)  +  iyi+i-yi+2) 

*  {yi-yi+2)  *  (yi-yi+i) 


(B.7x) 


A  =  In 


i^l  +  2~  ^k) 


(2/  -  ^k) 


(B.7y) 


B  = 


e{l  +  2,l;k) 


(B.7z) 


The  symbolic  algebra  software  package  DERIVE  was  of  great  help  in  the 
derivation  of  these  complicated  expressions. 
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